CN105549078B - The five dimension interpolation process methods and device of irregular seismic data - Google Patents

Abstract

Five the invention discloses a kind of irregular seismic data tie up interpolation process method and device, this method comprises: irregular seismic data is transformed to frequency domain and spatial domain；For the frequency slice of each seismic data for transforming to frequency domain and spatial domain, according to fast fourier transform algorithm, matrix and vector to seismic data carry out product calculation, the frequency-wavenumber domain data after seeking five dimension interpolation；Frequency-wavenumber domain data after five dimension interpolation are transformed into frequency domain and spatial domain, then transform to time-domain and spatial domain.Above-mentioned technical proposal realizes the quick processing to irregular seismic data, and the matrix of the entire treatment process overwhelming majority and vector product operation all utilize Fast Fourier Transform (FFT) fft algorithm to realize, improve the efficiency of seismic data process.

Description

The five dimension interpolation process methods and device of irregular seismic data

Technical field

The present invention relates to seismic exploration technique field, in particular to five dimension interpolation processing sides of a kind of irregular seismic data Method and device.

Background technique

Seismic prospecting is broadly divided into seismic data acquisition, processing and explains three phases.In the seismic data acquisition stage, nothing By being land exploration or seafari, influenced by various complicated factors, the spatial sampling of actual acquisition data can only be close It is even irregular like rule.For example, being explored for land, on the ground of the surface conditions such as city, highway, network of waterways complexity Area, can not layout rules observation system, it is necessary to irregularly acquired；It for seafari, is influenced by ocean current, practical electricity Cable survey line the phenomenon that there are featherings, equally it is difficult to ensure the spatial sampling of rule.In the seism processing stage, many numbers The earthquake of rule space sampling is required or at least had benefited from according to Processing Algorithm (such as migration algorithm and multiple removal algorithm) Data.In addition, carrying out in flakes or when time-lapse seismic data processing, the observation system ginseng of the seismic data of different batches acquisition Number is generally different, this also relates to the spatial sampling regularization problem of seismic data.Therefore, special interpolation algorithm is utilized It is very necessary for solving seismic data spatial sampling irregular problem.

Currently, most of seismic data interpolation algorithms are also limited only to three-dimensional interpolation, but the seismic data sheet of field acquisition It is the function of five dimension coordinates in matter, bidimensional is for determining shot point spatial position, and bidimensional is for determining geophone station spatial position, also One-dimensional to be used to determine the sampled point time, five dimension interpolation algorithms of development can be more fully using the seismic data of acquisition, and then obtains Obtain more preferable interpolation result.Compared with D interpolation algorithm, the primary problem that five dimension interpolation algorithms face is exactly huge data volume And calculation amount.The three-dimensional minimum weight norm interpolation algorithm of Liu and Sacchi (2004) are extended to five dimensions by Trad (2009), are calculated When method is realized, all matrixes and vector product operation may be by Fast Fourier Transform (FFT) (FFT) completion, hereby it is ensured that Computational efficiency, but algorithm require seismic data spatial sampling be it is regular, be not used to irregularly acquire seismic data, have Significant limitation.Jin (2010) is based on irregular spatial sampling it is assumed that proposing based on decaying minimum norm Fourier inverting Five dimension seismic data interpolation algorithms, algorithm using unequal interval Fast Fourier Transform (FFT) (NFFT) realize frequent matrix with to Product calculation is measured, improves computational efficiency to a certain extent, but since the computational efficiency of high dimensional data NFFT is far away from FFT, And NFFT itself is a kind of approximate algorithm, and the method for Jin is still to be improved in terms of computational efficiency.

Summary of the invention

Five the embodiment of the invention provides a kind of irregular seismic data tie up interpolation process method, to improve earthquake number According to the efficiency of processing, this method comprises:

Irregular seismic data is transformed into frequency domain and spatial domain；

For the frequency slice of each seismic data for transforming to frequency domain and spatial domain, according to Fast Fourier Transform (FFT) Algorithm, matrix and vector to seismic data carry out product calculation, the frequency-wavenumber domain data after seeking five dimension interpolation；

Frequency-wavenumber domain data after five dimension interpolation are transformed into frequency domain and spatial domain, then transform to time-domain and space Domain.

Five the embodiment of the invention provides a kind of irregular seismic data tie up interpolation processor, to improve earthquake number According to the efficiency of processing, which includes:

Seismic data conversion module, for irregular seismic data to be transformed to frequency domain and spatial domain；

Five dimension interpolation processing modules, cut for transforming to the frequency of seismic data of frequency domain and spatial domain for each Piece, according to fast fourier transform algorithm, matrix and vector to seismic data carry out product calculation, after seeking five dimension interpolation Frequency-wavenumber domain data；

Five dimension interpolated data conversion modules, for the frequency-wavenumber domain data after five dimension interpolation to be transformed to frequency domain and sky Between domain, then transform to time-domain and spatial domain.

Compared with prior art, technical solution provided in an embodiment of the present invention, by converting irregular seismic data To frequency domain and spatial domain；For the frequency slice of each seismic data for transforming to frequency domain and spatial domain, according to quick Fourier Transform Algorithm, matrix and vector to seismic data carry out product calculation, the frequency-wavenumber domain after seeking five dimension interpolation Data；Frequency-wavenumber domain data after five dimension interpolation are transformed into frequency domain and spatial domain, then transform to time-domain and spatial domain, The quick processing to irregular seismic data is realized, entire calculating process only relates to a small amount of NFFT operation, most squares Battle array may be by FFT realization with vector product operation, substantially increases computational efficiency, improves the effect of seismic data process Rate has important practical application value.

Detailed description of the invention

The drawings described herein are used to provide a further understanding of the present invention, constitutes part of this application, not Constitute limitation of the invention.In the accompanying drawings:

Fig. 1 is the flow diagram of five dimension interpolation process methods of irregular seismic data in the embodiment of the present invention；

Fig. 2 is the seismic data observation system schematic diagram in the embodiment of the present invention before interpolation；

Fig. 3 is the seismic data observation system schematic diagram of the same CMP face element in the embodiment of the present invention before interpolation；

Fig. 4 is the seismic data observation system of the same CMP face element in the embodiment of the present invention after interpolation corresponding with Fig. 3 System schematic diagram；

Fig. 5 is to indicate to be intended to using the operation time comparison diagram of different fast algorithms in the embodiment of the present invention；

Fig. 6 is the CMP trace gather seismic data schematic diagram in the embodiment of the present invention before interpolation corresponding with Fig. 3；

Fig. 7 is the CMP trace gather seismic data in the embodiment of the present invention after interpolation corresponding with Fig. 4；

Fig. 8 is the structural schematic diagram of five dimension interpolation processors of irregular seismic data in the embodiment of the present invention.

Specific embodiment

To make the objectives, technical solutions, and advantages of the present invention clearer, right below with reference to embodiment and attached drawing The present invention is described in further details.Here, exemplary embodiment and its explanation of the invention is used to explain the present invention, but simultaneously It is not as a limitation of the invention.

The invention proposes a kind of new iterative scheme, entire calculating process only relates to a small amount of NFFT operation, the overwhelming majority Matrix and vector product operation may be by FFT realization, substantially increase computational efficiency, have important practical application valence Value.1 to 8 pair of program describes in detail as follows with reference to the accompanying drawing.

Fig. 1 is the flow diagram of five dimension interpolation process methods of irregular seismic data in the embodiment of the present invention, such as Fig. 1 Shown, which includes the following steps:

Step 101: irregular seismic data is transformed into frequency domain and spatial domain；

Step 102: for the frequency slice of each seismic data for transforming to frequency domain and spatial domain, according to quick Fu In leaf transformation algorithm, matrix to seismic data and vector carry out product calculation, the frequency-wavenumber domain number after seeking five dimension interpolation According to；

Step 103: the frequency-wavenumber domain data after five dimension interpolation being transformed into frequency domain and spatial domain, then transform to the time Domain and spatial domain.

Compared with prior art, technical solution provided in an embodiment of the present invention, by converting irregular seismic data To frequency domain and spatial domain；For the frequency slice of each seismic data for transforming to frequency domain and spatial domain, according to quick Fourier Transform Algorithm, matrix and vector to seismic data carry out product calculation, the frequency-wavenumber domain after seeking five dimension interpolation Data；Frequency-wavenumber domain data after five dimension interpolation are transformed into frequency domain and spatial domain, then transform to time-domain and spatial domain, The quick processing to irregular seismic data is realized, entire calculating process only relates to a small amount of NFFT operation, most squares Battle array may be by FFT realization with vector product operation, substantially increases computational efficiency, improves the effect of seismic data process Rate has important practical application value.

In one embodiment, may include: before above-mentioned steps 101

Irregular seismic data is denoised, static correction and linear NMO are handled；

From the irregular seismic data extracted in treated irregular seismic data to five dimension interpolation processings；

When it is implemented, irregular seismic data is denoised, the processing of static correction and linear NMO, it may include: pair The seismic data of acquisition carries out preceding processing operations, including denoising, static correction, linear NMO etc..

When it is implemented, extracting the irregular earthquake number to five dimension interpolation processings from treated irregular seismic data According to may include: that extraction is desired carries out five seismic datas for tieing up interpolation

Above-mentioned t_itFor time sampling coordinate, the value range of subscript it is 0,1,2 ..., N_t- 1, i.e. time orientation has N_tIt is a Sampled point, if time sampling interval is Δ t, then t_it=it Δ t；

It is above-mentionedFor spatial sampling coordinate, the value range of subscript ix is 0,1,2 ..., N_x- 1, i.e. direction in space has N_xIt is a Sampled point,It is a point in space-time, whereinIt is that one kind typicallys represent form, according to actual needs, shot point x coordinate, shot point y-coordinate, detection can be respectively represented Point x coordinate, geophone station y-coordinate, can also respectively represent central point x coordinate, central point y-coordinate, geophone offset, azimuth.

Above-mentioned steps 101 may include: that the irregular seismic data to five dimension interpolation processings of extraction is transformed to frequency Domain and spatial domain.

When it is implemented, the irregular seismic data to five dimension interpolation processings of extraction is transformed to frequency domain and spatial domain May include:

It is rightIn every one-dimensional coordinate transformation is normalized, transformation results are denoted as x_ix=((x₀)_ix,(x₁)_ix,(x₂)_ix,(x₃)_ix), transformed seismic data is expressed as

Above-mentioned normalization transform method isWherein, η={ 0,1,2,3 } indicates not Same dimension, for any one-dimensional coordinateIts maximum value and minimum value are taken, is denoted as respectively vmax_ηAnd vmin_η, specificallyWithn_ηIt is according to reality Need and the coordinate points after the interpolation that sets, after interpolation the sampling interval of each space dimension be

Using Fast Fourier Transform (FFT) by seismic data d (x_ix,t_it) transform to frequency domain d (x_ix,ω_iω)。

Above-mentioned d (x_ix,ω_iω) and d (x_ix,t_it) meet following relational expression, Wherein,t_it=it Δ t, ω_iω=i ω Δ ω,I ω=0,1,2 ..., N_t-1}.Here Need to guarantee seismic data d (x_ix,t_it) meet sampling thheorem, if seismic data d (x_ix,t_it) useful signal energy just locate InFrequency band in, whereinThen need to meet condition i ω_max≤floor(N_t/2)。

Generate the discrete inverse Fourier transform matrix F (x of unequal interval_ix,k_ik), wherein ix={ 0,1,2 ..., N_xIt -1 } is square The line index of battle array, ik={ 0,1,2 ..., n₀n₁n₂n₃- 1 } it is indexed for matrix column.

It is above-mentionedWherein,It is in space-time One point,c_η=floor ((n_η- 1)/2),η=0,1,2,3, function Floor () indicates to be rounded downwards, k_ikx_ixIndicate four dimensional vector k_ikWith x_ixInner product.And ik with It keeps closing It is formulaWherein, n₀、n₁、n₂、n₃Definition see step (3).

In one embodiment, above-mentioned steps 102 may include:

To each frequency slice i ω=0,1,2 ..., i ω_max, it is handled 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 ω_maxIt is the upper frequency limit index chosen according to actual needs, selection principle is so that seismic data d (x_ix,t_it) useful signal energy just atFrequency band in；

w₀It include n for one₀n₁n₂n₃The column vector of a element, and w₀Each element be equal toIts In, w₀It is w_iωThe case where when middle subscript i ω=0；

F=F (x_ix,k_ik), it is a N_xRow, n₀n₁n₂n₃The matrix of column；Wherein, F (x_ix,k_ik) it is that unequal interval is discrete inverse Fourier transform matrix, ix={ 0,1,2 ..., N_xIt -1 } is the line index of matrix, ik={ 0,1,2 ..., n₀n₁n₂n₃It -1 } is square The column index of battle array,Wherein,It is in space-time One point,c_η=floor ((n_η- 1)/2),Function Floor () indicates to be rounded downwards, k_ikx_ixIndicate four dimensional vector k_ikWith x_ixInner product, ik withIt keeps closing It is formula

It include N for one_xThe column of a element Vector；

Function conj () indicates to seek each element of input vector complex conjugate, and exports a vector；Operator * It indicates that the corresponding element of two vectors is multiplied, and obtains a vector；Function sqrt () indicates each member to input vector Element extraction of square root, and export a vector；Function sum () indicates to sum to all elements of input vector, and exports a mark Amount；Operator^HThe conjugate transposition of matrix is sought in expression；

Function y=CG (b, w) is expressed as follows preconditioned conjugate gradient method, 1. sets constant, comprising: iteration ends are opposite Tolerance tol=10^-4, maximum number of iterations maxit=30, A=F^HF；2. calculating iteration initializaing variable, comprising:r₀=b, ρ_-1=1；3. calculating iteration ends absolute toleranceIt is right In i_cg={ 0,1,2 ..., maxit-1 }, executes following processing operation:

IfOr i_cg=maxit jumps out circulation, returnsOtherwise it continues cycling through:

Wherein,For scalar,Indicate scalarRespectively with vectorEach element multiplication, and A vector is returned,Indicate scalarRespectively with vectorEach element multiplication, and return a vector；

In one embodiment, in the function y=CG (b, w), each iteration needs to calculate a matrix and multiplies with vector Operation q=Ap is accumulated, wherein A=F^HF has Multilevel Block Toeplitz matrix structure, calculates q=Ap using following fast algorithm:

1. setting N₀=n₁n₂n₃、N₁=n₂n₃、N₂=n₃、N₃=1, m_η=2n_η- 1, η=0,1,2,3, it generates as follows Length is m_ηVectorWithWherein η=0,1,2,3:

2. setting vectorIt is a four-dimensional array Vectorization indicate, i_aFor the index of vector a,For the corresponding four-dimensional array indexing of vector a,Value range be 0,1,2,…,m_η- 1, η=0,1,2,3, hereinIndicate i_aIt isFunction, be all made of below Similar expression, as follows to the element assignment of a, initializing variable

ForExecute following circulation (1) to (12)

(1)

(2)

(3) forExecute following circulation (4) to (12)

(4)

(5)

(6) forExecute following circulation (7) to (12)

(7)

(8)

(9) forExecute following circulation (10) to (12)

(10)

(11)

(12)

Wherein,Indicate vector?A element, η=0,1,2,3,Indicate vector a's TheA element,The of representing matrix ARow,The element of column；

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

Wherein, vectorIt is a four-dimensional array Vectorization indicate,For vectorIndex,For the corresponding four-dimensional array indexing of vector a,Value range be 0,1,2,…,m_η- 1, η=0,1,2,3,Function FFT₄() expression does four-dimensional fast Fourier to four-dimensional array Transformation, and return to a four-dimensional array；

4. setting vectorIt is a four-dimensional array Vectorization indicate, i_pFor the index of vector p,For the corresponding four-dimensional array indexing of vector p,Value range be 0,1,2,…,n_η- 1, η=0,1,2,3；

If vectorIt is four dimensions The vectorization expression of group, i_p′For the index of vector p ',For the corresponding four-dimensional array indexing of vector p ',Value Range is 0,1,2 ..., m_η- 1, η=0,1,2,3；As follows to each element assignment of vector p ':

5. calculatingWherein, function FFT₄() and IFFT₄() respectively indicates to four dimensions Group does four-dimensional Fast Fourier Transform (FFT) and inverse transformation, and returns to a four-dimensional array；Operator * indicates two four-dimensional arrays Corresponding element is multiplied, and obtains a four-dimensional array；

VectorIt is a four-dimensional array Vectorization indicate, i_q′For the index of vector q ',For the corresponding four-dimensional array indexing of vector q ',Value model Enclosing is 0,1,2 ..., m_η- 1, η=0,1,2,3；

6. setting vectorIt is a four-dimensional array Vectorization indicate, i_qFor the index of vector q,For the corresponding four-dimensional array indexing of vector q,Value range be 0,1,2,…,n_η- 1, η=0,1,2,3, as follows to each element assignment of vector q,Wherein 0≤i_η≤n_η- 1, η=0,1,2,3.

When it is implemented, in entire method flow, above the 1.~3. step only need to calculate primary, operand substantially can be with Ignore, therefore, as described in 5. step, the calculation amount of matrix and vector product operation q=Ap mainly include primary four-dimensional quick Fu In leaf transformation and primary four-dimensional inverse fast Fourier transform.As shown in figure 5, the technical solution that the present inventor's embodiment provides improves The efficiency of irregular seismic data process.

In one embodiment, further includes: to y_iωExecute functionWherein, the value range of i ω is 0, 1,2,…,iω_max；

FunctionExpression is adjusted the element position of four-dimensional array y, and returns to a four-dimensional arrayIts In:

It is the vectorization of a four-dimensional array It indicates, i_yFor the index of vector y,For the corresponding four-dimensional array indexing of vector y,Value range be 0,1, 2,…,n_η- 1, η=0,1,2,3；

It is the vectorization of a four-dimensional array It indicates,For vectorIndex,For vectorCorresponding four-dimension array indexing,Value range be 0,1, 2,…,n_η- 1, η=0,1,2,3；

Whereinc_η=floor ((n_η-1)/ 2),Value range be 0,1,2 ..., n_η- 1, η=0,1,2,3, function mod (v₁,v₂) indicate v₁To v₂Seek 0~v₂Between -1 Remainder, function floor () indicate downwards be rounded.

In one embodiment, further includes: rightExecute functionWherein, the value range of i ω is 0,1,2,…,iω_max；

FunctionIt indicates to four-dimensional arrayFour-dimensional inverse fast Fourier transform is done, and returns to four dimensions Group s, wherein

It is the vectorization of a four-dimensional array It indicates,For vectorIndex,For vectorCorresponding four-dimension array indexing,Value range be 0,1, 2,…,n_η- 1, η=0,1,2,3；

It is the vectorization table of a four-dimensional array Show, i_sFor the index of vector s,For the corresponding four-dimensional array indexing of vector s, i_sValue range be 0,1,2 ..., n_η- 1, η=0,1,2,3.

In one embodiment, above-mentioned steps 103 may include:

Utilize s_iωFive dimension data of frequency domain after generating interpolationWherein,It is sat for the spatial sampling after interpolation Mark,Value range be 0,1,2 ..., n₀n₁n₂n₃The value upper limit of -1, i ω are by i ω_maxBecome N_t- 1, with meet it is subsequent into The demand of row inverse fast Fourier transform,

Wherein,i_ηValue range be 0,1, 2,…,n_η- 1, vmin_ηFor the minimum value of space dimension coordinate each after interpolation,Between sampling for space dimension each after interpolation Every η=0,1,2,3, ix and i₀、i₁、i₂、i₃Keep relational expression ix=((i₀·n₁+i₁)·n₂+i₂)·n₃+i₃； Generating mode is as follows, uses vectorIt representsFrequency coordinate index be i ω a part of data,

Then

Wherein, function conj () indicates to seek each element of input vector complex conjugate, and z is indicated by n₀n₁n₂n₃A 0 yuan The vector that element is constituted；

It will using inverse fast Fourier transformThe time-domain of transformation, five dimension datas after obtaining interpolation

Wherein,WithMeet following relational expression,Its In,t_itThe value range of=it Δ t, it are 0,1,2 ..., N_t- 1,Value range be 0,1,2 ..., n₀n₁n₂n₃- 1, N_tFor number of sampling points；Δ t is time sampling interval.

Based on the same inventive concept, a kind of five Wei Chazhichu of irregular seismic data are additionally provided in the embodiment of the present invention Device is managed, as described in the following examples.Due to irregular seismic data five dimension interpolation processor problems principle with not Five dimension interpolation process methods of regular seismic data are similar, therefore the implementation of five dimension interpolation processors of irregular seismic data It may refer to the implementation of five dimension interpolation process methods of irregular seismic data, overlaps will not be repeated.It is used below, The combination of the software and/or hardware of predetermined function may be implemented in term " unit " or " module ".Although following embodiment is retouched The device stated preferably realized with software, but the combined realization of hardware or software and hardware be also may and by structure Think.

Fig. 8 is the structural schematic diagram of five dimension interpolation processors of irregular seismic data in the embodiment of the present invention, such as Fig. 8 Shown, which includes:

Seismic data conversion module 10, for irregular seismic data to be transformed to frequency domain and spatial domain；

Five dimension interpolation processing modules 20, the frequency of the seismic data for transforming to frequency domain and spatial domain for each Slice, according to fast fourier transform algorithm, matrix and vector to seismic data carry out product calculation, after seeking five dimension interpolation Frequency-wavenumber domain data；

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

In an example, further includes:

Seismic data preprocessing module, for being denoised, at static correction and linear NMO to irregular seismic data Reason；

Seismic data abstraction module, for being extracted to five dimension interpolation processings not from treated irregular seismic data Regular seismic data；

The seismic data conversion module is specifically used for: the irregular seismic data to five dimension interpolation processings of extraction is become Change to frequency domain and spatial domain.

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

To each frequency slice i ω=0,1,2 ..., i ω_max, it is handled 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 ω_maxIt is the upper frequency limit index chosen according to actual needs, selection principle is so that seismic data d (x_ix,t_it) useful signal energy just atFrequency band in, wherein

w₀It include n for one₀n₁n₂n₃The column vector of a element, and w₀Each element be equal toIts In, w₀It is w_iωThe case where when middle subscript i ω=0；

F=F (x_ix,k_ik), it is a N_xRow, n₀n₁n₂n₃The matrix of column；Wherein, F (x_ix,k_ik) it is that unequal interval is discrete inverse Fourier transform matrix, ix are the line index of matrix, and the value range of ix is 0,1,2 ..., N_x- 1, ik are matrix column index, The value range of ik is 0,1,2 ..., n₀n₁n₂n₃- 1,Wherein,For imaginary unit,It is a point in space-time,c_η=floor ((n_η-1)/ 2),Value range be 0,1,2 ..., n_ηThe value range of -1, η are 0,1,2,3, and function floor () indicates to be rounded downwards, k_ikx_ixIndicate four dimensional vector k_ikWith x_ixInner product, ik with Keep relational expression

It include N for one_xThe column of a element Vector；

Function conj () indicates to seek each element of input vector complex conjugate, and exports a vector；Operator * It indicates that the corresponding element of two vectors is multiplied, and obtains a vector；Function sqrt () indicates each member to input vector Element extraction of square root, and export a vector；Function sum () indicates to sum to all elements of input vector, and exports a mark Amount；Operator^HThe conjugate transposition of matrix is sought in expression；

Function y=CG (b, w) is expressed as follows preconditioned conjugate gradient method, 1. sets constant, comprising: iteration ends are opposite Tolerance tol=10^-4, maximum number of iterations maxit=30, A=F^HF；2. calculating iteration initializaing variable, comprising:r₀=b, ρ_-1=1；3. calculating iteration ends absolute toleranceIt is right In i_cg={ 0,1,2 ..., maxit-1 }, executes following processing operation:

IfOr i_cg=maxit jumps out circulation, returnsOtherwise it continues cycling through:

Wherein,For scalar,Indicate scalarRespectively with vectorEach element multiplication, and A vector is returned,Indicate scalarRespectively with vectorEach element multiplication, and return a vector；

In an example, in the function y=CG (b, w), each iteration needs to calculate a matrix and vector product Operation q=Ap, wherein A=F^HF has Multilevel Block Toeplitz matrix structure, calculates q=Ap using following fast algorithm:

1. setting N₀=n₁n₂n₃、N₁=n₂n₃、N₂=n₃、N₃=1, m_η=2n_η- 1, η=0,1,2,3, it generates as follows Length is m_ηVectorWithWherein η=0,1,2,3:

2. setting vectorIt is a four-dimensional array Vectorization indicate, i_aFor the index of vector a,For the corresponding four-dimensional array indexing of vector a,Value range be 0,1,2,…,m_η- 1, η=0,1,2,3, hereinIndicate i_aIt isFunction, be all made of below Similar expression, as follows to the element assignment of a, initializing variable

ForExecute following circulation (1) to (12)

(1)

(2)

(3) forExecute following circulation (4) to (12)

(4)

(5)

(6) forExecute following circulation (7) to (12)

(7)

(8)

(9) forExecute following circulation (10) to (12)

(10)

(11)

(12)

Wherein,Indicate vector?A element, η=0,1,2,3,Indicate vector a's TheA element,The of representing matrix ARow,The element of column；

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

Wherein,Be a four-dimensional array to Quantization means,For vectorIndex,For the corresponding four-dimensional array indexing of vector a,Value range be 0,1, 2,…,m_η- 1, η=0,1,2,3,Function FFT₄() expression does four-dimensional Fast Fourier Transform (FFT) to four-dimensional array, And return to a four-dimensional array；

4. setting vectorIt is a four-dimensional array Vectorization indicate, i_pFor the index of vector p,For the corresponding four-dimensional array indexing of vector p,Value range be 0,1,2,…,n_η- 1, η=0,1,2,3；

If vectorIt is four dimensions The vectorization expression of group, i_p′For the index of vector p ',For the corresponding four-dimensional array indexing of vector p ',Value Range is 0,1,2 ..., m_η- 1, η=0,1,2,3；As follows to each element assignment of vector p ':

5. calculatingWherein, function FFT₄() and IFFT₄() respectively indicates to four dimensions Group does four-dimensional Fast Fourier Transform (FFT) and inverse transformation, and returns to a four-dimensional array；Operator * indicates two four-dimensional arrays Corresponding element is multiplied, and obtains a four-dimensional array；

Wherein, vectorIt is one four The vectorization expression of dimension group, i_q′For the index of vector q ',For the corresponding four-dimensional array indexing of vector q ','s Value range is 0,1,2 ..., m_η- 1, η=0,1,2,3；

6. setting vectorIt is a four-dimensional array Vectorization indicate, i_qFor the index of vector q,For the corresponding four-dimensional array indexing of vector q,Value range be 0,1,2,…,n_η- 1, η=0,1,2,3, as follows to each element assignment of vector q,Wherein 0≤i_η≤n_η- 1, η=0,1,2,3.

In an example, further includes: to y_iωExecute functionWherein, the value range of i ω is 0,1, 2,…,iω_max；

FunctionExpression is adjusted the element position of four-dimensional array y, and returns to a four-dimensional arrayIts In:

It is the vectorization of a four-dimensional array It indicates, i_yFor the index of vector y,For the corresponding four-dimensional array indexing of vector y,Value range be 0,1, 2,…,n_η- 1, η=0,1,2,3；

It is the vectorization of a four-dimensional array It indicates,For vectorIndex,For vectorCorresponding four-dimension array indexing,Value range be 0,1, 2,…,n_η- 1, η=0,1,2,3；

Whereinc_η=floor ((n_η-1)/ 2),Value range be 0,1,2 ..., n_η- 1, η=0,1,2,3, function mod (v₁,v₂) indicate v₁To v₂Seek 0~v₂Between -1 Remainder, function floor () indicate downwards be rounded.

In an example, further includes: rightExecute functionWherein, the value range of i ω is 0, 1,2,…,iω_max；

FunctionIt indicates to four-dimensional arrayFour-dimensional inverse fast Fourier transform is done, and returns to four dimensions Group s, wherein

It is the vectorization of a four-dimensional array It indicates,For vectorIndex,For vectorCorresponding four-dimension array indexing,Value range be 0,1, 2,…,n_η- 1, η=0,1,2,3；

It is the vectorization table of a four-dimensional array Show, i_sFor the index of vector s,For the corresponding four-dimensional array indexing of vector s, i_sValue range be 0,1,2 ..., n_η- 1, η=0,1,2,3.

In an example, the five dimensions interpolated data conversion module is specifically used for:

Utilize s_iωFive dimension data of frequency domain after generating interpolationWherein,It is sat for the spatial sampling after interpolation Mark,Value range be 0,1,2 ..., n₀n₁n₂n₃The value upper limit of -1, i ω are by i ω_maxBecome N_t- 1, with meet it is subsequent into The demand of row inverse fast Fourier transform,

Wherein,i_ηValue range be 0,1, 2,…,n_η- 1, vmin_ηFor the minimum value of space dimension coordinate each after interpolation,Between sampling for space dimension each after interpolation Every η=0,1,2,3, ix and i₀、i₁、i₂、i₃Keep relational expression ix=((i₀·n₁+i₁)·n₂+i₂)·n₃+i₃； Generating mode is as follows, uses vectorIt representsFrequency coordinate index be i ω a part of data,

Then

Wherein, function conj () indicates to seek each element of input vector complex conjugate, and z is indicated by n₀n₁n₂n₃A 0 yuan The vector that element is constituted；

It will using inverse fast Fourier transformThe time-domain of transformation, five dimension datas after obtaining interpolation

Wherein,WithMeet following relational expression,Its In,t_itThe value range of=it Δ t, it are 0,1,2 ..., N_t- 1,Value range be 0,1,2 ..., n₀n₁n₂n₃- 1, N_tFor number of sampling points；Δ t is time sampling interval.

It is illustrated again with example below, in order to understanding how to implement the present invention.

The five dimension interpolation process methods that the embodiment of the present invention provides irregular seismic data mainly include the following steps:

(1), preceding processing operations, including denoising, static correction, linear NMO etc. are carried out to the seismic data of acquisition；

(2), the seismic data for wanting to carry out five dimension interpolation is extractedFig. 2 is interpolation provided in an embodiment of the present invention Preceding seismic data layout chart.

T described in step (2)_itFor time sampling coordinate, subscript it={ 0,1,2 ..., N_t- 1 }, i.e., time orientation has N_t A 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 }, i.e., direction in space has N_x A sampled point,It is a point in space-time, whereinIt is that one kind typicallys represent form, according to this example needs, respectively represents central point x coordinate, central point y-coordinate, big gun inspection Away from azimuth.In this example, N_x=7151.

(3), rightIn every one-dimensional coordinate transformation, transformation results are normalized It is denoted as x_ix=((x₀)_ix,(x₁)_ix,(x₂)_ix,(x₃)_ix), transformed seismic data is expressed as

Normalization transform method described in step (3) isWherein, η=0,1, 2,3 } different dimensions is indicated, for any one-dimensional coordinateIts maximum value and minimum value are taken, It is denoted as vmax respectively_ηAnd vmin_η, specificallyWithn_ηIt is Coordinate points after the interpolation set according to actual needs, after interpolation the sampling interval of each space dimension beIn 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,

(4), using Fast Fourier Transform (FFT) by seismic data d (x_ix,t_it) transform to frequency domain d (x_ix,ω_iω)。

D (x described in step (4)_ix,ω_iω) and d (x_ix,t_it) meet following relational expression,Wherein,t_it=it Δ t, ω_iω=i ω Δ ω,I ω=0,1,2 ..., N_t-1}.It needs exist for guaranteeing seismic data d (x_ix,t_it) meet sampling thheorem, if ground Shake data d (x_ix,t_it) useful signal energy just atFrequency band in, whereinThen need Meet condition i ω_max≤floor(N_t/2).In this example, Δ ω ≈ 2.618, i ω_max=175.

(5), the discrete inverse Fourier transform matrix F (x of unequal interval is generated_ix,k_ik), wherein ix={ 0,1,2 ..., N_x-1} For the line index of matrix, ik={ 0,1,2 ..., n₀n₁n₂n₃- 1 } it is indexed for matrix column.In this example, n₀n₁n₂n₃=141570.

Step (5) is describedWherein,It is four-dimensional empty Between in a point,c_η=floor ((n_η- 1)/2),η=0,1,2,3, letter Number floor () indicates to be rounded downwards, k_ikx_ixIndicate four dimensional vector k_ikWith x_ixInner product.And ik withIt keeps Relational expressionWherein, 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, following circulation is executed, and store y_iω。

I ω described in step (6)_maxIt is the upper frequency limit index chosen according to actual needs, selection principle is so that earthquake Data d (x_ix,t_it) useful signal energy just atFrequency band in, wherein

W described in step (6)₀It include n for one₀n₁n₂n₃The column vector of a element, and w₀Each element be equal toWherein, w₀It is w_iωThe case where when middle subscript i ω=0.

F=F (x described in step (6)_ix,k_ik), it is a N_xRow, n₀n₁n₂n₃The matrix of column.

Described in step (6)Include for one N_xThe column vector of a element.

Function conj () described in step (6) indicates to seek each element of input vector complex conjugate, and exports one Vector；Operator * indicates that the corresponding element of two vectors is multiplied, and obtains a vector；Function sqrt () is indicated to input Each element of vector extracts square root, and exports a vector；Function sum () indicates to sum to all elements of input vector, And export a scalar；Operator^HThe conjugate transposition of matrix is sought in expression, similarly hereinafter.

Function y=CG (b, w) described in step (6) is expressed as follows preconditioned conjugate gradient method, 1. sets constant.Including Iteration ends relative tolerance tol=10^-4, maximum number of iterations maxit=30, A=F^HF；2. calculating iteration initializaing variable.Includingr₀=b, ρ_-1=1；3. calculating iteration ends absolute toleranceFor i_cg={ 0,1,2 ..., maxit } executes following circulation

IfOr i_cg=maxit jumps out circulation, returnsOtherwise it continues cycling through:

Wherein,For scalar,Indicate scalarRespectively with vectorEach element multiplication, and return A vector is returned,Indicate scalarRespectively with vectorEach element multiplication, and return a vector；

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

1. setting N₀=n₁n₂n₃、N₁=n₂n₃、N₂=n₃、N₃=1, m_η=2n_η- 1, η=0,1,2,3, it generates as follows Length is m_ηVectorWithWherein η=0,1,2,3:

2. setting vectorIt is a four-dimensional array Vectorization indicate, i_aFor the index of vector a,For the corresponding four-dimensional array indexing of vector a,Value range be 0,1,2,…,m_η- 1, η=0,1,2,3, hereinIndicate i_aIt isFunction, be all made of below Similar expression.As follows to the element assignment of a.Initializing variable

ForExecute following circulation (1) to (12)

(1)

(2)

(3) forExecute following circulation (4) to (12)

(4)

(5)

(6) forExecute following circulation (7) to (12)

(7)

(8)

(9) forExecute following circulation (10) to (12)

(10)

(11)

(12)

Wherein,Indicate vector?A element, η=0,1,2,3,Indicate vector a's TheA element,The of representing matrix ARow,The element of column.

3. doing four-dimensional Fast Fourier Transform (FFT), and return vector to four-dimensional array represented by vector aWherein, vector Wherein,It is the vectorization table of a four-dimensional array Show,For vectorIndex,For the corresponding four-dimensional array indexing of vector a,Value range be 0,1,2 ..., m_η- 1, η=0,1,2,3.Function FFT₄() expression does four-dimensional Fast Fourier Transform (FFT) to four-dimensional array, and returns Return a four-dimensional array.

4. setting vectorIt is a four-dimensional array Vectorization indicate, i_pFor the index of vector p,For the corresponding four-dimensional array indexing of vector p,Value range be 0,1,2,…,n_η- 1, η=0,1,2,3；

If vectorIt is four dimensions The vectorization expression of group, i_p′For the index of vector p ',For the corresponding four-dimensional array indexing of vector p ',Value Range is 0,1,2 ..., m_η- 1, η=0,1,2,3.As follows to each element assignment of vector p '

5. calculatingWherein, function FFT₄() and IFFT₄() respectively indicates to four dimensions Group does four-dimensional Fast Fourier Transform (FFT) and inverse transformation, and returns to a four-dimensional array；Operator * indicates two four-dimensional arrays Corresponding element is multiplied, and obtains a four-dimensional array；Wherein, vectorIt is the vectorization table of a four-dimensional array Show, i_q′For the index of vector q ',For the corresponding four-dimensional array indexing of vector q ',Value range be 0,1, 2,…,m_η- 1, η=0,1,2,3.

6. setting vectorIt is a four-dimensional array Vectorization indicate, i_qFor the index of vector q,For the corresponding four-dimensional array indexing of vector q,Value range be 0,1,2,…,n_η- 1, η=0,1,2,3.As follows to each element assignment of vector q,Wherein 0≤i_η≤n_η- 1, η=0,1,2,3.

7. in entire method flow, above the 1.~3. step only need to calculate primary, operand can be ignored substantially, because This, as described in 5. step, the calculation amount of matrix and vector product operation q=Ap mainly include primary four-dimensional Fast Fourier Transform (FFT) With primary four-dimensional inverse fast Fourier transform.

(7), to y_iωExecute functionAnd it storesWherein the value range of i ω is 0,1,2 ..., i ω_max。

Function described in step (7)Expression is adjusted the element position of four-dimensional array y, and returns to one A four-dimension arrayWherein,

It is the vectorization of a four-dimensional array It indicates, i_yFor the index of vector y,For the corresponding four-dimensional array indexing of vector y,Value range be 0,1, 2,…,n_η- 1, η=0,1,2,3；

It is the vectorization of a four-dimensional array It indicates,For vectorIndex,For vectorCorresponding four-dimension array indexing,Value range be 0,1, 2,…,n_η- 1, η=0,1,2,3.

Specific method of adjustment is as follows,Wherein c_η=floor ((n_η- 1)/2),Value range be 0,1,2 ..., n_η- 1, η=0,1,2,3, function mod (v₁,v₂) indicate v₁To v₂Seek 0~v₂Remainder between -1, function floor () indicate to be rounded downwards.

(8), rightExecute functionAnd store s_iω, wherein the value range of i ω is 0,1,2 ..., iω_max。

Function described in step (8)It indicates to four-dimensional arrayFour-dimensional inverse fast Fourier transform is done, and Return to a four-dimensional array s.Wherein,

It is the vectorization of a four-dimensional array It indicates,For vectorIndex,For vectorCorresponding four-dimension array indexing,Value range be 0,1, 2,…,n_η- 1, η=0,1,2,3；

It is the vectorization table of a four-dimensional array Show, i_sFor the index of vector s,For the corresponding four-dimensional array indexing of vector s, i_sValue range be 0,1,2 ..., n_η- 1, η=0,1,2,3.

(9), s is utilized_iωFive dimension data of frequency domain after generating interpolationWherein,It is adopted for the space after interpolation Sample coordinate,Value range be 0,1,2 ..., n₀n₁n₂n₃The value upper limit of -1, i ω are by i ω_maxBecome N_t- 1, after meeting The continuous demand for carrying out inverse fast Fourier transform.

Described in step (9)

Wherein,i_ηValue range be 0,1, 2,…,n_η- 1, vmin_ηFor the minimum value of space dimension coordinate each after interpolation,Step (3) are shown in definition.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, states for convenience, uses vectorIt represents A part of data, specifically,Then

Wherein, function conj () expression seeks complex conjugate to each element of input vector.Z is indicated by n₀n₁n₂n₃A 0 yuan The vector that element is constituted.

It (10), will using inverse fast Fourier transformThe time-domain of transformation, five dimension datas after obtaining interpolation

Described in step (10)WithMeet following relational expression,Wherein,t_itThe value range of=it Δ t, it are 0,1, 2,…,N_t- 1,Value range be 0,1,2 ..., n₀n₁n₂n₃- 1, N_tFor number of sampling points；Δ t is time sampling interval.

Fig. 2 is the seismic data observation system schematic diagram in the embodiment of the present invention before interpolation；Fig. 3 is in the embodiment of the present invention The seismic data observation system schematic diagram of the same CMP face element before interpolation；Fig. 4 is corresponding with Fig. 3 in the embodiment of the present invention The seismic data observation system schematic diagram of the same CMP face element after interpolation；Fig. 5 is in the embodiment of the present invention using different fast The operation time comparison diagram of the short-cut counting method indicates to be intended to；Fig. 6 is the CMP trace gather in the embodiment of the present invention before interpolation corresponding with Fig. 3 Seismic data schematic diagram；Fig. 7 is the CMP trace gather seismic data in the embodiment of the present invention after interpolation corresponding with Fig. 4.

Compared with prior art, by shown in Fig. 2 to Fig. 7 and the introduction of above-described embodiment it is found that the embodiment of the present invention Five dimension interpolation fast algorithms of the irregular acquisition seismic data provided, core are to propose a kind of new iterative scheme, significantly Improve computational efficiency.Entire calculating process only relates to a small amount of NFFT operation, and most matrixes and vector product operation are all It can use FFT realization, substantially increase computational efficiency, there is important practical application value.

It should be understood by those skilled in the art that, the embodiment of the present invention can provide as method, system or computer program Product.Therefore, complete hardware embodiment, complete software embodiment or reality combining software and hardware aspects can be used in the present invention Apply the form of example.Moreover, it wherein includes the computer of computer usable program code that the present invention, which can be used in one or more, The computer program implemented in usable storage medium (including but not limited to magnetic disk storage, CD-ROM, optical memory etc.) produces The form of product.

The present invention be referring to according to the method for the embodiment of the present invention, the process of equipment (system) and computer program product Figure and/or block diagram describe.It should be understood that every one stream in flowchart and/or the block diagram can be realized by computer program instructions The combination of process and/or box in journey and/or box and flowchart and/or the block diagram.It can provide these computer programs Instruct the processor of general purpose computer, special purpose computer, Embedded Processor or other programmable data processing devices to produce A raw machine, so that being generated by the instruction that computer or the processor of other programmable data processing devices execute for real The device for the function of being specified in present one or more flows of the flowchart and/or one or more blocks of the block diagram.

These computer program instructions, which may also be stored in, is able to guide computer or other programmable data processing devices with spy Determine in the computer-readable memory that mode works, so that it includes referring to that instruction stored in the computer readable memory, which generates, Enable the manufacture of device, the command device realize in one box of one or more flows of the flowchart and/or block diagram or The function of being specified in multiple boxes.

These computer program instructions also can be loaded onto a computer or other programmable data processing device, so that counting Series of operation steps are executed on calculation machine or other programmable devices to generate computer implemented processing, thus in computer or The instruction executed on other programmable devices is provided for realizing in one or more flows of the flowchart and/or block diagram one The step of function of being specified in a box or multiple boxes.

The foregoing is only a preferred embodiment of the present invention, is not intended to restrict the invention, for the skill of this field For art personnel, the embodiment of the present invention can have various modifications and variations.All within the spirits and principles of the present invention, made Any modification, equivalent substitution, improvement and etc. should all be included in the protection scope of the present invention.

Claims (10)

1. a kind of five dimension interpolation process methods of irregular seismic data characterized by comprising

Irregular seismic data is transformed into frequency domain and spatial domain；

For the frequency slice of each seismic data for transforming to frequency domain and spatial domain, calculated according to Fast Fourier Transform (FFT) Method, matrix and vector to the seismic data for transforming to frequency domain and spatial domain carry out product calculation, after seeking five dimension interpolation Frequency-wavenumber domain data；

Frequency-wavenumber domain data after five dimension interpolation are transformed into frequency domain and spatial domain, then transform to time-domain and spatial domain；

For the frequency slice of each seismic data for transforming to frequency domain and spatial domain, calculated according to Fast Fourier Transform (FFT) Method, matrix and vector to the seismic data for transforming to frequency domain and spatial domain carry out product calculation, after seeking five dimension interpolation Frequency-wavenumber domain data, comprising:

To each frequency slice i ω=0,1,2 ..., i ω_max, it is handled 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 ω_maxIt is the upper frequency limit index chosen according to actual needs, selection principle is so that seismic data d (x_ix,t_it) Useful signal energy just atFrequency band in；

w₀It include n for one₀n₁n₂n₃The column vector of a element, and w₀Each element be equal toWherein, w₀ It is w_iωThe case where when middle subscript i ω=0；

F=F (x_ix,k_ik), it is a N_xRow, n₀n₁n₂n₃The matrix of column；Wherein, F (x_ix,k_ik) it is in discrete inverse Fu of unequal interval Leaf transformation matrix, ix are the line index of matrix, and the value range of ix is 0,1,2 ..., N_x- 1, ik are matrix column index, ik's Value range is 0,1,2 ..., n₀n₁n₂n₃- 1,Wherein,For imaginary unit,It is a point in space-time,c_η=floor ((n_η-1)/ 2),Value range be 0,1,2 ..., n_ηThe value range of -1, η are 0,1,2,3, and function floor () indicates to be rounded downwards, k_ikx_ixIndicate four dimensional vector k_ikWith x_ixInner product；

It include N for one_xThe column vector of a element；

Function conj () indicates to seek each element of input vector complex conjugate, and exports a vector；Operator * is indicated The corresponding element of two vectors is multiplied, and obtains a vector；Each element of input vector is opened in function sqrt () expression Square, and export a vector；Function sum () indicates to sum to all elements of input vector, and exports a scalar；Fortune Operator H indicates to seek the conjugate transposition of 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 number of iterations maxit=30, A=F^HF；2. calculating iteration initializaing variable, comprising:r₀=b, ρ_-1=1；3. calculating iteration ends absolute toleranceIt is right In i_cg={ 0,1,2 ..., maxit-1 }, executes following processing operation:

IfOr i_cg=maxit jumps out circulation, returnsOtherwise it continues cycling through:

Wherein,For scalar,Indicate scalarRespectively with vectorEach element multiplication, and return one A vector,Indicate scalarRespectively with vectorEach element multiplication, and return a vector；

x_ixFor the spatial sampling coordinate of transformed seismic data, t_itFor the time sampling coordinate of seismic data；

d(x_ix,ω_iω) it is d (x_ix,t_it) to transform to the frequency after frequency domain be ω_iωSeismic data, ω_iωSubscript i ω be The discrete index of frequency coordinate, d_iωFor frequency slice vector；b_iωFor frequency slice vector d_iωUnequal interval Fourier transformation； y_iωFor the output of preconditioned conjugate gradient method；u_iωFor y_iωTwo norm distances square；w_iωFor pre-conditional conjugate gradient calculation The input weight of method；

In the function y=CG (b, w), each iteration needs to calculate a matrix and vector product operation q=Ap, wherein A= F^HF has Multilevel Block Toeplitz matrix structure, calculates q=Ap using following fast algorithm:

1. setting N₀=n₁n₂n₃、N₁=n₂n₃、N₂=n₃、N₃=1, m_η=2n_η- 1, η=0,1,2,3, length is generated as follows For m_ηVectorWithWherein η=0,1,2,3:

2. setting vectorBe a four-dimensional array to Quantization means, i_aFor the index of vector a,For the corresponding four-dimensional array indexing of vector a,Value range be 0,1, 2,…,m_η- 1, η=0,1,2,3, hereinIndicate i_aIt isFunction, be all made of class below Like expression, as follows to the element assignment of a, initializing variable

ForExecute following circulation (1) to (12)

(3) forExecute following circulation (4) to (12)

(6) forExecute following circulation (7) to (12)

(9) forExecute following circulation (10) to (12)

Wherein,Indicate vector?A element, η=0,1,2,3,Indicate the of vector aA element,The of representing matrix ARow,The element of column；

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

Wherein, vectorBe a four-dimensional array to Quantization means,For vectorIndex,For the corresponding four-dimensional array indexing of vector a,Value range be 0,1, 2,…,m_η- 1, η=0,1,2,3, function FFT₄() indicates to do four-dimensional array four-dimensional Fast Fourier Transform (FFT), and returns to one Four-dimensional array；

4. setting vectorBe a four-dimensional array to Quantization means, i_pFor the index of vector p,For the corresponding four-dimensional array indexing of vector p,Value range be 0,1, 2,…,n_η- 1, η=0,1,2,3；

If vectorIt is a four-dimensional array Vectorization expression, i_p′For the index of vector p ',For the corresponding four-dimensional array indexing of vector p ',Value range It is 0,1,2 ..., m_η- 1, η=0,1,2,3；As follows to each element assignment of vector p ':

5. calculatingWherein, function FFT₄() and IFFT₄(), which respectively indicates, does four to four-dimensional array Fast Fourier Transform (FFT) and inverse transformation are tieed up, and returns to a four-dimensional array；Operator * indicates the corresponding element of two four-dimensional arrays Element is multiplied, and obtains a four-dimensional array；

VectorBe a four-dimensional array to Quantization means, i_q′For the index of vector q ',For the corresponding four-dimensional array indexing of vector q ',Value range be 0,1,2,…,m_η- 1, η=0,1,2,3；

6. setting vectorBe a four-dimensional array to Quantization means, i_qFor the index of vector q,For the corresponding four-dimensional array indexing of vector q,Value range be 0,1, 2,…,n_η- 1, η=0,1,2,3, as follows to each element assignment of vector q,Wherein i_ηValue range be 0,1,2 ..., n_η- 1。

2. five dimension interpolation process methods of irregular seismic data as described in claim 1, which is characterized in that described to advise Before then seismic data transforms to frequency domain and spatial domain, comprising:

Irregular seismic data is denoised, the processing of static correction and linear NMO；

From the irregular seismic data extracted in treated irregular seismic data to five dimension interpolation processings；

It is described that irregular seismic data is transformed into frequency domain and spatial domain, comprising: not to five dimension interpolation processings by extraction Regular seismic data transforms to frequency domain and spatial domain.

3. five dimension interpolation process methods of irregular seismic data as described in claim 1, which is characterized in that further include: it is right y_iωExecute functionWherein, the value range of i ω is 0,1,2 ..., i ω_max；

FunctionExpression is adjusted the element position of four-dimensional array y, and returns to a four-dimensional arrayWherein:

It is the vectorization expression of a four-dimensional array, i_yFor the index of vector y,For the corresponding four-dimensional array indexing of vector y,Value range be 0,1,2 ..., n_η- 1, η=0,1,2,3；

It is the vectorization expression of a four-dimensional array,For vectorIndex,For vectorCorresponding four-dimension array indexing,Value range be 0,1,2 ..., n_η- 1, η=0,1,2,3；

Whereinc_η=floor ((n_η- 1)/2), Value range be 0,1,2 ..., n_η- 1, η=0,1,2,3, function mod (v₁,v₂) indicate v₁To v₂Seek 0~v₂It is remaining between -1 Number, function floor () indicate to be rounded downwards.

4. five dimension interpolation process methods of irregular seismic data as claimed in claim 3, which is characterized in that further include: it is rightExecute functionWherein, the value range of i ω is 0,1,2 ..., i ω_max；

FunctionIt indicates to four-dimensional arrayFour-dimensional inverse fast Fourier transform is done, and returns to a four-dimensional array s, Wherein,

It is the vectorization table of a four-dimensional array Show,For vectorIndex,For vectorCorresponding four-dimension array indexing,Value range be 0,1,2 ..., n_η- 1, η=0,1,2,3；

It is the vectorization expression of a four-dimensional array, i_s For the index of vector s,For the corresponding four-dimensional array indexing of vector s, i_sValue range be 0,1,2 ..., n_η- 1, η =0,1,2,3.

5. five dimension interpolation process methods of irregular seismic data as claimed in claim 4, which is characterized in that after interpolation Frequency-wavenumber domain data transform to frequency and spatial domain, then transform to time and spatial domain, complete to irregular seismic data Five dimension interpolation processings, comprising:

Utilize s_iωFive dimension data of frequency domain after generating interpolationWherein,For the spatial sampling coordinate after interpolation,Value range be 0,1,2 ..., n₀n₁n₂n₃The value upper limit of -1, i ω are by i ω_maxBecome N_t- 1, to meet subsequent carry out fastly The demand of fast inverse Fourier transform；

Wherein,i_ηValue range be 0,1,2 ..., n_η- 1, vmin_ηFor the minimum value of space dimension coordinate each after interpolation,For the sampling interval of space dimension each after interpolation, η= 0,1,2,3；Generating mode is as follows, uses vectorIt representsFrequency coordinate index be one of i ω Divided data,Then

Wherein, function conj () indicates to seek each element of input vector complex conjugate, and z is indicated by n₀n₁n₂n₃A 0 element structure At vector；

It will using inverse fast Fourier transformThe time-domain of transformation, five dimension datas after obtaining interpolation

Wherein,WithMeet following relational expression,Wherein,t_itThe value range of=it Δ t, it are 0,1,2 ..., N_t- 1,Value range be 0,1,2 ..., n₀n₁n₂n₃- 1, N_tFor number of sampling points；Δ t is time sampling interval.

6. a kind of five dimension interpolation processors of irregular seismic data characterized by comprising

Seismic data conversion module, for irregular seismic data to be transformed to frequency domain and spatial domain；

Five dimension interpolation processing modules, the frequency slice of the seismic data for transforming to frequency domain and spatial domain for each, According to fast fourier transform algorithm, matrix and vector to the seismic data for transforming to frequency domain and spatial domain carry out product fortune It calculates, the frequency-wavenumber domain data after seeking five dimension interpolation；

Five dimension interpolated data conversion modules, for the frequency-wavenumber domain data after five dimension interpolation to be transformed to frequency domain and space Domain, then transform to time-domain and spatial domain；

The five dimensions interpolation processing module is specifically used for:

To each frequency slice i ω=0,1,2 ..., i ω_max, it is handled 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 ω_maxIt is the upper frequency limit index chosen according to actual needs, selection principle is so that seismic data d (x_ix,t_it) Useful signal energy just atFrequency band in, wherein

w₀It include n for one₀n₁n₂n₃The column vector of a element, and w₀Each element be equal toWherein, w₀ It is w_iωThe case where when middle subscript i ω=0；

F=F (x_ix,k_ik), it is a N_xRow, n₀n₁n₂n₃The matrix of column；Wherein, F (x_ix,k_ik) it is in discrete inverse Fu of unequal interval Leaf transformation matrix, ix are the line index of matrix, and the value range of ix is 0,1,2 ..., N_x- 1, ik are matrix column index, ik's Value range is 0,1,2 ..., n₀n₁n₂n₃- 1,Wherein,For imaginary unit,It is a point in space-time,c_η=floor ((n_η-1)/ 2),Value range be 0,1,2 ..., n_ηThe value range of -1, η are 0,1,2,3, and function floor () indicates to be rounded downwards, k_ikx_ixIndicate four dimensional vector k_ikWith x_ixInner product；

It include N for one_xThe column vector of a element；

Function conj () indicates to seek each element of input vector complex conjugate, and exports a vector；Operator * is indicated The corresponding element of two vectors is multiplied, and obtains a vector；Each element of input vector is opened in function sqrt () expression Square, and export a vector；Function sum () indicates to sum to all elements of input vector, and exports a scalar；Fortune Operator H indicates to seek the conjugate transposition of 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 number of iterations maxit=30, A=F^HF；2. calculating iteration initializaing variable, comprising:r₀=b, ρ_-1=1；3. calculating iteration ends absolute toleranceIt is right In i_cg={ 0,1,2 ..., maxit-1 }, executes following processing operation:

IfOr i_cg=maxit jumps out circulation, returnsOtherwise it continues cycling through:

Wherein,For scalar,Indicate scalarRespectively with vectorEach element multiplication, and return one A vector,Indicate scalarRespectively with vectorEach element multiplication, and return a vector；

x_ixFor the spatial sampling coordinate of transformed seismic data, t_itFor the time sampling coordinate of seismic data；

d(x_ix,ω_iω) it is d (x_ix,t_it) to transform to the frequency after frequency domain be ω_iωSeismic data, ω_iωSubscript i ω be The discrete index of frequency coordinate, d_iωFor frequency slice vector；b_iωFor frequency slice vector d_iωUnequal interval Fourier transformation； y_iωFor the output of preconditioned conjugate gradient method；u_iωFor y_iωTwo norm distances square；w_iωFor pre-conditional conjugate gradient calculation The input weight of method；

In the function y=CG (b, w), each iteration needs to calculate a matrix and vector product operation q=Ap, wherein A= F^HF has Multilevel Block Toeplitz matrix structure, calculates q=Ap using following fast algorithm:

1. setting N₀=n₁n₂n₃、N₁=n₂n₃、N₂=n₃、N₃=1, m_η=2n_η- 1, η=0,1,2,3, length is generated as follows For m_ηVectorWithWherein η=0,1,2,3:

2. setting vectorBe a four-dimensional array to Quantization means, i_aFor the index of vector a,For the corresponding four-dimensional array indexing of vector a,Value range be 0,1, 2,…,m_η- 1, η=0,1,2,3, hereinIndicate i_aIt isFunction, be all made of class below Like expression, as follows to the element assignment of a, initializing variable

ForExecute following circulation (1) to (12)

(3) forExecute following circulation (4) to (12)

(6) forExecute following circulation (7) to (12)

(9) forExecute following circulation (10) to (12)

Wherein,Indicate vector?A element, η=0,1,2,3,Indicate the of vector aA element,The of representing matrix ARow,The element of column；

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

Wherein,It is the vectorization of a four-dimensional array It indicates,For vectorIndex,For the corresponding four-dimensional array indexing of vector a,Value range be 0,1, 2,…,m_η- 1, η=0,1,2,3,Function FFT₄() expression does four-dimensional Fast Fourier Transform (FFT) to four-dimensional array, And return to a four-dimensional array；

4. setting vectorBe a four-dimensional array to Quantization means, i_pFor the index of vector p,For the corresponding four-dimensional array indexing of vector p,Value range be 0,1, 2,…,n_η- 1, η=0,1,2,3；

If vectorIt is a four-dimensional array Vectorization expression, i_p′For the index of vector p ',For the corresponding four-dimensional array indexing of vector p ',Value range It is 0,1,2 ..., m_η- 1, η=0,1,2,3；As follows to each element assignment of vector p ':

5. calculatingWherein, function FFT₄() and IFFT₄(), which respectively indicates, does four to four-dimensional array Fast Fourier Transform (FFT) and inverse transformation are tieed up, and returns to a four-dimensional array；Operator * indicates the corresponding element of two four-dimensional arrays Element is multiplied, and obtains a four-dimensional array；

VectorBe a four-dimensional array to Quantization means, i_q′For the index of vector q ',For the corresponding four-dimensional array indexing of vector q ',Value range be 0,1,2,…,m_η- 1, η=0,1,2,3；

6. setting vectorBe a four-dimensional array to Quantization means, i_qFor the index of vector q,For the corresponding four-dimensional array indexing of vector q,Value range be 0,1, 2,…,n_η- 1, η=0,1,2,3, as follows to each element assignment of vector q,Wherein i_ηValue range be 0,1,2 ..., n_η- 1。

7. five dimension interpolation processors of irregular seismic data as claimed in claim 6, which is characterized in that further include:

Seismic data preprocessing module, for being denoised to irregular seismic data, static correction and linear NMO processing；

Seismic data abstraction module, for extracting from treated irregular seismic data to the irregular of five dimension interpolation processings Seismic data；

The seismic data conversion module is specifically used for: the irregular seismic data to five dimension interpolation processings of extraction is transformed to Frequency domain and spatial domain.

8. five dimension interpolation processors of irregular seismic data as claimed in claim 6, which is characterized in that further include: it is right y_iωExecute functionWherein, the value range of i ω is 0,1,2 ..., i ω_max；

FunctionExpression is adjusted the element position of four-dimensional array y, and returns to a four-dimensional arrayWherein:

It is the vectorization expression of a four-dimensional array, i_yFor the index of vector y,For the corresponding four-dimensional array indexing of vector y,Value range be 0,1,2 ..., n_η- 1, η=0,1,2,3；

It is the vectorization expression of a four-dimensional array,For vectorIndex,For vectorCorresponding four-dimension array indexing,Value range be 0,1,2 ..., n_η- 1, η=0,1,2,3；

Whereinc_η=floor ((n_η- 1)/2), Value range be 0,1,2 ..., n_η- 1, η=0,1,2,3, function mod (v₁,v₂) indicate v₁To v₂Seek 0~v₂It is remaining between -1 Number, function floor () indicate to be rounded downwards.

9. five dimension interpolation processors of irregular seismic data as claimed in claim 8, which is characterized in that further include:

It is rightExecute functionWherein, the value range of i ω is 0,1,2 ..., i ω_max；

FunctionIt indicates to four-dimensional arrayFour-dimensional inverse fast Fourier transform is done, and returns to a four-dimensional array s, Wherein,

It is the vectorization expression of a four-dimensional array,For vectorIndex,For vectorCorresponding four-dimension array indexing,Value range be 0,1,2 ..., n_η- 1, η=0,1,2,3；

It is the vectorization expression of a four-dimensional array, i_s For the index of vector s,For the corresponding four-dimensional array indexing of vector s, i_sValue range be 0,1,2 ..., n_η- 1, η =0,1,2,3.

10. five dimension interpolation processors of irregular seismic data as claimed in claim 9, which is characterized in that five dimension Interpolated data conversion module is specifically used for:

Utilize s_iωFive dimension data of frequency domain after generating interpolationWherein,For the spatial sampling coordinate after interpolation,Value range be 0,1,2 ..., n₀n₁n₂n₃The value upper limit of -1, i ω are by i ω_maxBecome N_t- 1, to meet subsequent carry out fastly The demand of fast inverse Fourier transform；

Wherein,i_ηValue range be 0,1,2 ..., n_η- 1, vmin_ηFor the minimum value of space dimension coordinate each after interpolation,For the sampling interval of space dimension each after interpolation, η= 0,1,2,3；Generating mode is as follows, uses vectorIt representsFrequency coordinate index be i ω a part Data,

Then

Wherein, function conj () indicates to seek each element of input vector complex conjugate, and z is indicated by n₀n₁n₂n₃A 0 element structure At vector；

It will using inverse fast Fourier transformThe time-domain of transformation, five dimension datas after obtaining interpolation

Wherein,WithMeet following relational expression,Wherein,t_itThe value range of=it Δ t, it are 0,1,2 ..., N_t- 1,Value range be 0,1,2 ..., n₀n₁n₂n₃- 1, N_tFor number of sampling points；Δ t is time sampling interval.