CN105549078B - The five dimension interpolation process methods and device of irregular seismic data - Google Patents
The five dimension interpolation process methods and device of irregular seismic data Download PDFInfo
- Publication number
- CN105549078B CN105549078B CN201511029487.7A CN201511029487A CN105549078B CN 105549078 B CN105549078 B CN 105549078B CN 201511029487 A CN201511029487 A CN 201511029487A CN 105549078 B CN105549078 B CN 105549078B
- Authority
- CN
- China
- Prior art keywords
- vector
- dimensional array
- seismic data
- value range
- frequency
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 230000001788 irregular Effects 0.000 title claims abstract description 62
- 238000000034 method Methods 0.000 title claims abstract description 58
- 230000008569 process Effects 0.000 title claims abstract description 25
- 239000013598 vector Substances 0.000 claims abstract description 276
- 239000011159 matrix material Substances 0.000 claims abstract description 56
- 238000012545 processing Methods 0.000 claims abstract description 40
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 29
- 238000004364 calculation method Methods 0.000 claims abstract description 16
- 230000001131 transforming effect Effects 0.000 claims abstract description 13
- 238000005070 sampling Methods 0.000 claims description 39
- 230000009466 transformation Effects 0.000 claims description 21
- 238000006243 chemical reaction Methods 0.000 claims description 10
- 238000013139 quantization Methods 0.000 claims description 10
- 238000000605 extraction Methods 0.000 claims description 8
- 238000002939 conjugate gradient method Methods 0.000 claims description 7
- 238000012937 correction Methods 0.000 claims description 7
- 230000003068 static effect Effects 0.000 claims description 7
- 241001269238 Data Species 0.000 claims description 6
- 238000003491 array Methods 0.000 claims description 5
- 230000001351 cycling effect Effects 0.000 claims description 5
- 230000017105 transposition Effects 0.000 claims description 5
- 238000007781 pre-processing Methods 0.000 claims description 2
- 230000006870 function Effects 0.000 description 40
- 238000010586 diagram Methods 0.000 description 20
- 238000004590 computer program Methods 0.000 description 7
- 230000000694 effects Effects 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 238000010606 normalization Methods 0.000 description 2
- 241000208340 Araliaceae Species 0.000 description 1
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 1
- 235000003140 Panax quinquefolius Nutrition 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000005611 electricity Effects 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 235000008434 ginseng Nutrition 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000005012 migration Effects 0.000 description 1
- 238000013508 migration Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/40—Transforming data representation
- G01V2210/48—Other transforms
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
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
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 titFor time sampling coordinate, the value range of subscript it is 0,1,2 ..., Nt- 1, i.e. time orientation has NtIt is a
Sampled point, if time sampling interval is Δ t, then tit=it Δ t;
It is above-mentionedFor spatial sampling coordinate, the value range of subscript ix is 0,1,2 ..., Nx- 1, i.e. direction in space has NxIt 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
xix=((x0)ix,(x1)ix,(x2)ix,(x3)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 (xix,tit) transform to frequency domain d (xix,ωiω)。
Above-mentioned d (xix,ωiω) and d (xix,tit) meet following relational expression,
Wherein,tit=it Δ t, ωiω=i ω Δ ω,I ω=0,1,2 ..., Nt-1}.Here
Need to guarantee seismic data d (xix,tit) meet sampling thheorem, if seismic data d (xix,tit) useful signal energy just locate
InFrequency band in, whereinThen need to meet condition i ωmax≤floor(Nt/2)。
Generate the discrete inverse Fourier transform matrix F (x of unequal intervalix,kik), wherein ix={ 0,1,2 ..., NxIt -1 } is square
The line index of battle array, ik={ 0,1,2 ..., n0n1n2n3- 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, kikxixIndicate four dimensional vector kikWith xixInner product.And ik with It keeps closing
It is formulaWherein, n0、n1、n2、n3Definition 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:
biω=FHdiω;
yiω=CG (biω,wiω);
uiω=conj (yiω)·*yiω;
wiω+1=(sqrt (uiω/sum(uiω))+wiω)/2;
Wherein, i ωmaxIt is the upper frequency limit index chosen according to actual needs, selection principle is so that seismic data d
(xix,tit) useful signal energy just atFrequency band in;
w0It include n for one0n1n2n3The column vector of a element, and w0Each element be equal toIts
In, w0It is wiωThe case where when middle subscript i ω=0;
F=F (xix,kik), it is a NxRow, n0n1n2n3The matrix of column;Wherein, F (xix,kik) it is that unequal interval is discrete inverse
Fourier transform matrix, ix={ 0,1,2 ..., NxIt -1 } is the line index of matrix, ik={ 0,1,2 ..., n0n1n2n3It -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, kikxixIndicate four dimensional vector kikWith xixInner product, ik withIt keeps closing
It is formula
It include N for onexThe 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;OperatorHThe 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=FHF;2. calculating iteration initializaing variable, comprising:r0=b, ρ-1=1;3. calculating iteration ends absolute toleranceIt is right
In icg={ 0,1,2 ..., maxit-1 }, executes following processing operation:
IfOr icg=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=FHF has Multilevel Block Toeplitz matrix structure, calculates q=Ap using following fast algorithm:
1. setting N0=n1n2n3、N1=n2n3、N2=n3、N3=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, iaFor 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 iaIt 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 FFT4() 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, ipFor 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, ip′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 FFT4() and IFFT4() 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, iq′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, iqFor 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 yiω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, iyFor 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 (v1,v2) indicate v1To v2Seek 0~v2Between -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, isFor the index of vector s,For the corresponding four-dimensional array indexing of vector s, isValue range be 0,1,2 ...,
nη- 1, η=0,1,2,3.
In one embodiment, above-mentioned steps 103 may include:
Utilize siω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 ..., n0n1n2n3The value upper limit of -1, i ω are by i ωmaxBecome Nt- 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 i0、i1、i2、i3Keep relational expression ix=((i0·n1+i1)·n2+i2)·n3+i3;
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 n0n1n2n3A 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,titThe value range of=it Δ t, it are 0,1,2 ..., Nt- 1,Value range be 0,1,2 ...,
n0n1n2n3- 1, NtFor 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:
biω=FHdiω;
yiω=CG (biω,wiω);
uiω=conj (yiω)·*yiω;
wiω+1=(sqrt (uiω/sum(uiω))+wiω)/2;
Wherein, i ωmaxIt is the upper frequency limit index chosen according to actual needs, selection principle is so that seismic data d
(xix,tit) useful signal energy just atFrequency band in, wherein
w0It include n for one0n1n2n3The column vector of a element, and w0Each element be equal toIts
In, w0It is wiωThe case where when middle subscript i ω=0;
F=F (xix,kik), it is a NxRow, n0n1n2n3The matrix of column;Wherein, F (xix,kik) 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 ..., Nx- 1, ik are matrix column index,
The value range of ik is 0,1,2 ..., n0n1n2n3- 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,
kikxixIndicate four dimensional vector kikWith xixInner product, ik with Keep relational expression
It include N for onexThe 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;OperatorHThe 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=FHF;2. calculating iteration initializaing variable, comprising:r0=b, ρ-1=1;3. calculating iteration ends absolute toleranceIt is right
In icg={ 0,1,2 ..., maxit-1 }, executes following processing operation:
IfOr icg=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=FHF has Multilevel Block Toeplitz matrix structure, calculates q=Ap using following fast algorithm:
1. setting N0=n1n2n3、N1=n2n3、N2=n3、N3=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, iaFor 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 iaIt 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 FFT4() 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, ipFor 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, ip′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 FFT4() and IFFT4() 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, iq′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, iqFor 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 yiω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, iyFor 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 (v1,v2) indicate v1To v2Seek 0~v2Between -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, isFor the index of vector s,For the corresponding four-dimensional array indexing of vector s, isValue 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 siω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 ..., n0n1n2n3The value upper limit of -1, i ω are by i ωmaxBecome Nt- 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 i0、i1、i2、i3Keep relational expression ix=((i0·n1+i1)·n2+i2)·n3+i3;
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 n0n1n2n3A 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,titThe value range of=it Δ t, it are 0,1,2 ..., Nt- 1,Value range be 0,1,2 ...,
n0n1n2n3- 1, NtFor 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 ..., Nt- 1 }, i.e., time orientation has Nt
A sampled point, if time sampling interval is Δ t, then tit=it Δ t.In this example, Nt=1200, Δ t=0.002 seconds;
Described in step (2)For spatial sampling coordinate, subscript ix={ 0,1,2 ..., Nx- 1 }, i.e., direction in space has Nx
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, Nx=7151.
(3), rightIn every one-dimensional coordinate transformation, transformation results are normalized
It is denoted as xix=((x0)ix,(x1)ix,(x2)ix,(x3)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, vmin0=253, vmax0=263, vmin1=660, vmax1=670, vmin2=-
160, vmax2=160, vmin3=0, vmax3=170, n0=11, n1=11, n2=65, n3=18,
(4), using Fast Fourier Transform (FFT) by seismic data d (xix,tit) transform to frequency domain d (xix,ωiω)。
D (x described in step (4)ix,ωiω) and d (xix,tit) meet following relational expression,Wherein,tit=it Δ t, ωiω=i ω Δ ω,I ω=0,1,2 ..., Nt-1}.It needs exist for guaranteeing seismic data d (xix,tit) meet sampling thheorem, if ground
Shake data d (xix,tit) useful signal energy just atFrequency band in, whereinThen need
Meet condition i ωmax≤floor(Nt/2).In this example, Δ ω ≈ 2.618, i ωmax=175.
(5), the discrete inverse Fourier transform matrix F (x of unequal interval is generatedix,kik), wherein ix={ 0,1,2 ..., Nx-1}
For the line index of matrix, ik={ 0,1,2 ..., n0n1n2n3- 1 } it is indexed for matrix column.In this example, n0n1n2n3=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, kikxixIndicate four dimensional vector kikWith xixInner product.And ik withIt keeps
Relational expressionWherein, n0、n1、n2、n3Definition see step (3).In this example, c0=5, c1
=5, c2=32, c3=8.
(6), to each frequency slice i ω=0,1,2 ..., i ωmax, following circulation is executed, and store yiω。
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 (xix,tit) useful signal energy just atFrequency band in, wherein
W described in step (6)0It include n for one0n1n2n3The column vector of a element, and w0Each element be equal toWherein, w0It is wiωThe case where when middle subscript i ω=0.
F=F (x described in step (6)ix,kik), it is a NxRow, n0n1n2n3The matrix of column.
Described in step (6)Include for one
NxThe 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;OperatorHThe 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=FHF;2. calculating iteration initializaing variable.Includingr0=b, ρ-1=1;3. calculating iteration ends absolute toleranceFor
icg={ 0,1,2 ..., maxit } executes following circulation
IfOr icg=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=FHF has Multilevel Block Toeplitz matrix structure, can calculate q=Ap using following fast algorithm.
1. setting N0=n1n2n3、N1=n2n3、N2=n3、N3=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, iaFor 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 iaIt 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 FFT4() 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, ipFor 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, ip′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 FFT4() and IFFT4() 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, iq′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, iqFor 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 yiω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, iyFor 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 (v1,v2) indicate
v1To v2Seek 0~v2Remainder between -1, function floor () indicate to be rounded downwards.
(8), rightExecute functionAnd store siω, 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, isFor the index of vector s,For the corresponding four-dimensional array indexing of vector s, isValue range be 0,1,2 ...,
nη- 1, η=0,1,2,3.
(9), s is utilizediω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 ..., n0n1n2n3The value upper limit of -1, i ω are by i ωmaxBecome Nt- 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 i0、i1、
i2、i3Keep relational expression ix=((i0·n1+i1)·n2+i2)·n3+i3。
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 n0n1n2n3A 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,titThe value range of=it Δ t, it are 0,1,
2,…,Nt- 1,Value range be 0,1,2 ..., n0n1n2n3- 1, NtFor 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:
biω=FHdiω;
yiω=CG (biω,wiω);
uiω=conj (yiω)·*yiω;
wiω+1=(sqrt (uiω/sum(uiω))+wiω)/2;
Wherein, i ωmaxIt is the upper frequency limit index chosen according to actual needs, selection principle is so that seismic data d (xix,tit)
Useful signal energy just atFrequency band in;
w0It include n for one0n1n2n3The column vector of a element, and w0Each element be equal toWherein, w0
It is wiωThe case where when middle subscript i ω=0;
F=F (xix,kik), it is a NxRow, n0n1n2n3The matrix of column;Wherein, F (xix,kik) 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 ..., Nx- 1, ik are matrix column index, ik's
Value range is 0,1,2 ..., n0n1n2n3- 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,
kikxixIndicate four dimensional vector kikWith xixInner product;
It include N for onexThe 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=FHF;2. calculating iteration initializaing variable, comprising:r0=b, ρ-1=1;3. calculating iteration ends absolute toleranceIt is right
In icg={ 0,1,2 ..., maxit-1 }, executes following processing operation:
IfOr icg=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;
xixFor the spatial sampling coordinate of transformed seismic data, titFor the time sampling coordinate of seismic data;
d(xix,ωiω) it is d (xix,tit) to transform to the frequency after frequency domain be ωiωSeismic data, ωiωSubscript i ω be
The discrete index of frequency coordinate, diωFor frequency slice vector;biωFor frequency slice vector diωUnequal interval Fourier transformation;
yiωFor the output of preconditioned conjugate gradient method;uiωFor yiωTwo norm distances square;wiω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=
FHF has Multilevel Block Toeplitz matrix structure, calculates q=Ap using following fast algorithm:
1. setting N0=n1n2n3、N1=n2n3、N2=n3、N3=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, iaFor 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 iaIt 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 FFT4() 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, ipFor 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, ip′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 FFT4() and IFFT4(), 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, iq′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, iqFor 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
yiω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,
iyFor 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 (v1,v2) indicate v1To v2Seek 0~v2It 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, is
For the index of vector s,For the corresponding four-dimensional array indexing of vector s, isValue 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 siωFive dimension data of frequency domain after generating interpolationWherein,For the spatial sampling coordinate after interpolation,Value range be 0,1,2 ..., n0n1n2n3The value upper limit of -1, i ω are by i ωmaxBecome Nt- 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 n0n1n2n3A 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,titThe value range of=it Δ t, it are 0,1,2 ..., Nt- 1,Value range be 0,1,2 ...,
n0n1n2n3- 1, NtFor 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:
biω=FHdiω;
yiω=CG (biω,wiω);
uiω=conj (yiω)·*yiω;
wiω+1=(sqrt (uiω/sum(uiω))+wiω)/2;
Wherein, i ωmaxIt is the upper frequency limit index chosen according to actual needs, selection principle is so that seismic data d (xix,tit)
Useful signal energy just atFrequency band in, wherein
w0It include n for one0n1n2n3The column vector of a element, and w0Each element be equal toWherein, w0
It is wiωThe case where when middle subscript i ω=0;
F=F (xix,kik), it is a NxRow, n0n1n2n3The matrix of column;Wherein, F (xix,kik) 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 ..., Nx- 1, ik are matrix column index, ik's
Value range is 0,1,2 ..., n0n1n2n3- 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,
kikxixIndicate four dimensional vector kikWith xixInner product;
It include N for onexThe 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=FHF;2. calculating iteration initializaing variable, comprising:r0=b, ρ-1=1;3. calculating iteration ends absolute toleranceIt is right
In icg={ 0,1,2 ..., maxit-1 }, executes following processing operation:
IfOr icg=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;
xixFor the spatial sampling coordinate of transformed seismic data, titFor the time sampling coordinate of seismic data;
d(xix,ωiω) it is d (xix,tit) to transform to the frequency after frequency domain be ωiωSeismic data, ωiωSubscript i ω be
The discrete index of frequency coordinate, diωFor frequency slice vector;biωFor frequency slice vector diωUnequal interval Fourier transformation;
yiωFor the output of preconditioned conjugate gradient method;uiωFor yiωTwo norm distances square;wiω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=
FHF has Multilevel Block Toeplitz matrix structure, calculates q=Ap using following fast algorithm:
1. setting N0=n1n2n3、N1=n2n3、N2=n3、N3=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, iaFor 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 iaIt 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 FFT4() 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, ipFor 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, ip′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 FFT4() and IFFT4(), 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, iq′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, iqFor 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
yiω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,
iyFor 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 (v1,v2) indicate v1To v2Seek 0~v2It 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, is
For the index of vector s,For the corresponding four-dimensional array indexing of vector s, isValue 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 siωFive dimension data of frequency domain after generating interpolationWherein,For the spatial sampling coordinate after interpolation,Value range be 0,1,2 ..., n0n1n2n3The value upper limit of -1, i ω are by i ωmaxBecome Nt- 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 n0n1n2n3A 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,titThe value range of=it Δ t, it are 0,1,2 ..., Nt- 1,Value range be 0,1,2 ...,
n0n1n2n3- 1, NtFor number of sampling points;Δ t is time sampling interval.
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201511029487.7A CN105549078B (en) | 2015-12-31 | 2015-12-31 | The five dimension interpolation process methods and device of irregular seismic data |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201511029487.7A CN105549078B (en) | 2015-12-31 | 2015-12-31 | The five dimension interpolation process methods and device of irregular seismic data |
Publications (2)
Publication Number | Publication Date |
---|---|
CN105549078A CN105549078A (en) | 2016-05-04 |
CN105549078B true CN105549078B (en) | 2019-06-11 |
Family
ID=55828378
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201511029487.7A Active CN105549078B (en) | 2015-12-31 | 2015-12-31 | The five dimension interpolation process methods and device of irregular seismic data |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN105549078B (en) |
Families Citing this family (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106646612B (en) * | 2016-12-20 | 2018-11-30 | 中国地质大学(北京) | Reconstruction of seismic data method based on matrix contraction |
CN107807389B (en) * | 2017-09-18 | 2019-07-09 | 中国石油天然气股份有限公司 | The seismic data encryption method and device of anti-alias |
CN107703539B (en) * | 2017-09-18 | 2019-05-07 | 中国石油天然气股份有限公司 | The seismic data interpolation method and device of anti-alias |
CN107561588B (en) * | 2017-09-19 | 2019-07-09 | 中国石油天然气股份有限公司 | A kind of seismic data noise drawing method and device |
CN108345034B (en) * | 2018-02-06 | 2021-08-03 | 北京中科海讯数字科技股份有限公司 | Seismic data regularization method |
CN112666608B (en) * | 2019-10-15 | 2024-06-25 | 中国石油天然气股份有限公司 | Five-dimensional interpolation method and system for steep dip seismic signals |
CN111929725B (en) * | 2020-07-29 | 2021-07-02 | 中国石油大学(北京) | Seismic data interpolation method, device and equipment |
Family Cites Families (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20120215453A1 (en) * | 2011-02-22 | 2012-08-23 | Cggveritas Services Sa | Device and method for multi-dimensional coherency driven denoising data |
US9158016B2 (en) * | 2012-04-30 | 2015-10-13 | Conocophillips Company | Multi-dimensional data reconstruction constrained by a regularly interpolated model |
CN104459770B (en) * | 2013-09-24 | 2017-06-16 | 中国石油化工股份有限公司 | A kind of method for regularizing high-dimensional seismic data |
-
2015
- 2015-12-31 CN CN201511029487.7A patent/CN105549078B/en active Active
Also Published As
Publication number | Publication date |
---|---|
CN105549078A (en) | 2016-05-04 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN105549078B (en) | The five dimension interpolation process methods and device of irregular seismic data | |
CN107153216B (en) | Determine the method, apparatus and computer storage medium of the Poynting vector of seismic wave field | |
CN106371138B (en) | Reconstruction of seismic data method and apparatus | |
Honkonen et al. | Predicting global ground geoelectric field with coupled geospace and three‐dimensional geomagnetic induction models | |
Jia et al. | A fast rank-reduction algorithm for three-dimensional seismic data interpolation | |
Liu et al. | A new kind of optimal second-order symplectic scheme for seismic wave simulations | |
Zhang et al. | 2-D seismic data reconstruction via truncated nuclear norm regularization | |
CN115292973B (en) | Arbitrarily sampled space wave number domain three-dimensional magnetic field numerical simulation method and system | |
CN114239268B (en) | Method for acquiring cross-interface radiation field of underwater double-electric-dipole array based on Romberg | |
Gao et al. | A fast rank reduction method for the reconstruction of 5D seismic volumes | |
Gao et al. | Fast least-squares reverse time migration via a superposition of Kronecker products | |
CN107356971B (en) | Seismic data rule method and device | |
Ceniceros et al. | A fast, robust, and non-stiff immersed boundary method | |
CN111291316A (en) | Multi-scale resistivity inversion method and system based on wavelet transformation | |
CN107144881B (en) | The treating method and apparatus of seismic data | |
Da Silva et al. | Applications of low-rank compressed seismic data to full-waveform inversion and extended image volumes | |
CN104166795B (en) | A kind of multiple sine wave frequency estimating methods based on many observation vector rarefaction representations | |
CN107703539B (en) | The seismic data interpolation method and device of anti-alias | |
Mastryukov | Optimal finite difference schemes for a wave equation | |
Chai et al. | Modeling multisource multifrequency acoustic wavefields by a multiscale Fourier feature physics-informed neural network with adaptive activation functions | |
Kumar et al. | Time-jittered marine acquisition: A rank-minimization approach for 5D source separation | |
Goodrick | Image reconstruction in radio astronomy with non-coplanar synthesis arrays | |
CN107807389B (en) | The seismic data encryption method and device of anti-alias | |
Yu et al. | Simultaneous off-the-grid regularization and reconstruction for 3D seismic data by a new combined sampling operator | |
Zhang et al. | Vandermonde constrained CANDECOMP/PARAFAC tensor decomposition for high-dimensional seismic data reconstruction |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |