Summary of the invention
To existing situation, the objective of the invention is: a kind of three repair and maintenance width of cloth prestack time migration methods are provided, through improve migration aperture, provide the weight coefficient of protecting the width of cloth, in imaging process compacting skew noise, realize efficient, high s/n ratio, guarantor's width of cloth migration imaging.This method can obtain to protect common reflection point (CRP) the road collection of the width of cloth, utilizes prestack information to carry out technology such as oil gas and fluid detection thereby serve prestack inversion etc. better.
The present invention takes following technical scheme:
(1) writes down the seismic reflection signals of artificial epicenter excitation with many towing cables or many surveys line, record on the tape;
(2) read seismic data from tape, the pre-stack seismic data is carried out conventional compacting noise treatment, to several common midpoints position; Extract CMP gather; These road collection are done conventional NMO correction (NMO) velocity pick, the gained result is done laterally on average, as initial, horizontal migration velocity field uniformly; The pre-stack seismic data is pressed offset distance ordering, is divided into groups, with not on the same group seismic data be placed on the various computing node of cluster computer;
(3) according to the inclination maximum of the offset distance of seismic data, azimuthal variation range and underground structure and initial, horizontal migration velocity field uniformly, the migration aperture table in the three-D migration aperture that becomes during structure;
(4) confirm corresponding compacting skew noise alleviation coefficient table according to the migration aperture table;
(5) weight coefficient that adopts the deconvolution image-forming condition based on one way wave equation and depth shift the to obtain amplitude calculating of squinting;
(6) at each computing node; To whole seismic trace circulations,, confirm the initial imaging time of each discrete point calculations of offset in effective imaging region and the zone by the migration aperture table to each seismic trace; Begin to calculate the skew amplitude of each imaging point from initial imaging time; To the imaging point in the marginarium that is in the three-D migration aperture, look into the attenuation coefficient table and pick up attenuation coefficient and multiply by the skew amplitude and obtain final amplitude, amplitude is added in the array of depositing migration result on the corresponding offset distance;
(7) collect the calculations of offset result of each computing node, to each horizontal level of imaging region the calculations of offset result by the time degree of depth, offset distance ordering, form the CRP gather of each horizontal position;
(8) whether straight according to the lineups in the CRP gather, in migration process, confirm the migration velocity field;
(9) utilize the migration velocity field of confirming in the migration process to do calculations of offset; Whole CRP gathers to forming are done RNMO; The corresponding numerical value that obviously stretching and noise part occur is changed to zero (i.e. excision), with the numerical value stack of different offset distances, formation skew stacked section;
(10) be the profile image of underground reflective construct through the software for display stacked section data-switching that will squint, profile image will be indicated form, breakpoint, turn-off size and the interface of the underground structure reflection strength to seismic event.
The inclination maximum of described offset distance, azimuthal variation range and underground structure according to seismic data, and initial laterally uniformly the migration aperture table in migration velocity field, the three-D migration aperture that becomes when making up be achieved in that with one group of data (x
i, y
i, 2T
i) (i=1 k) describes the three-D migration aperture, and wherein k is total discrete counting in effective imaging region (imaging space is in the projection of surface level), (x
i, y
i) be that two horizontal coordinates of certain discrete point and seismic trace central point are poor in the imaging region, and T
iIt then is the initial imaging time (during outward journey) of this calculations of offset; Making the shot point coordinate of seismic trace is (x
s, y
s), the geophone station coordinate is (x
r, y
r), the migration aperture that becomes during introducing is exactly in the calculations of offset of this seismic trace, is (x to horizontal coordinate in the imaging region
i+ 0.5 (x
s+ x
r), y
i+ 0.5 (y
s+ y
r)) point (i=1, k), from T=T
iBegin to carry out calculations of offset.
By two inclination maximums of structure along survey line and vertical line direction, θ
xAnd θ
y, make up following inequation group:
In the formula
Be the position angle of seismic trace, satisfy the minimum angles α and the β of last two formulas, guarantee that exactly inclination maximum is θ
xAnd θ
yStructure when being able to correctly form images, along two maximum imaging angles of azimuth direction and vertical orientations angular direction.
To different time degree of depth T; It is that long axis direction, vertical orientations angular direction are short-axis direction, the elliptical area of center on the central point of seismic trace that imaging region can be regarded as with the azimuth direction, by α, β and initial, laterally migration velocity field
can be tried to achieve its long and short semiaxis and is respectively uniformly
To be defined as
be half offset distance of seismic trace to parameter lambda in the formula, when time depth T is outward journey.
To time depth T and T+ Δ T (Δ T is the sampling interval that forms images on the time depth direction); Calculate corresponding long and short semiaxis respectively; Two ellipses that make this two group leader, minor semi-axis form have common central point (0; 0), these two oval bands that form be exactly initial imaging time be the imaging region of T; With the discrete point coordinate x=m Δ x that drops in this band, it is right that y=n Δ y and T constitute data, and wherein Δ x and Δ y are the horizontal coordinate samplings of migration imaging.To time depth sampled point circulation (greater than this time depth, migration aperture will not change in time), to the whole (x that obtain less than the fixed time degree of depth; Y; The size of 2T) pressing x and y respectively sorts, and adds up the k that always counts, the three-D migration aperture that becomes in the time of just can obtaining.
The position angle and the offset distance of seismic trace have determined migration aperture; Need position angle and offset distance variation range according to whole seismic traces; (be the minimizing memory space according to angle and offset distance sampling interval; Desirable bigger spacing), calculates the three-D migration aperture that changes with position angle and offset distance, make up the migration aperture table; The migration aperture table is a three-dimensional matrice, the corresponding offset distance of dimension, and another dimension counterparty parallactic angle, the third dimension are three coordinate data (x in three-D migration aperture
i, y
i, 2T
i) and n
iTotally 4 numbers increase arrangement, n here with i
iIt is one group of integer of corresponding attenuation coefficient; First element of this one dimension is the k that always counts of this migration aperture.That in fact, deposit in this matrix is coordinate (x
i, y
i, 2T
i) corresponding discrete point sequence number, so the migration aperture table is an INTEGER MATRICES.
Azimuthal sampling interval is not to the actual angle value, but to azimuthal sine value.When picking up migration aperture, examine sine value, offset distance and position angle sine value are rounded by sampling interval, directly in table, pick up in the coordinate difference computer azimuth angle of vertical side line direction by the offset distance and the big gun of seismic trace.
Described compacting skew noise alleviation coefficient table according to the definite correspondence of migration aperture table is achieved in that the corresponding attenuation coefficient array of each migration aperture is obtained by following steps:
1) the definition three-dimensional matrice (x, y, T), according to (the x that obtains from migration aperture
i, y
i) point initial imaging time T
b, design factor
And deposit it in correspondence (x
i, y
i, T) in the position, n is counting of comprising of the marginarium of artificial definition in the formula, Δ T is the sampling interval (during outward journey) of time depth direction imaging, to all k dot cycles in the migration aperture;
2) to each time depth T, by computes a and b:
Each parameter identical when calculating migration aperture in the formula; By the horizontal direction imaging sampling interval Δ y of maximum, calculating
and then calculate
wherein,
is the position angle of seismic trace.(x to f>1
i, y
i, T), if corresponding (x
i, y
i, T) stored c in the position
i, do not do operation, otherwise calculate
With c
2Deposit in the corresponding position;
3) (each element T) circulates, and will in preceding two steps, not have the some assignment 1.0 of assignment for x, y to three-dimensional matrice; It is level and smooth that matrix element is space Gauss, only keeps the element less than 1.0 after level and smooth, forms the attenuation coefficient array by these elements.Every group of (x of migration aperture array
i, y
i, 2T
i) corresponding sequence c in the attenuation coefficient array
j(j=1, n
i), n wherein
iBe matrix (x, y, T) in (x
i, y
i) locate, along T
bPlay element number, c less than 1.0
jBe these less than 1.0 element value, as discussing in the relevant migration aperture table, n
iTo leave in the migration aperture table.
The attenuation coefficient table is corresponding with the migration aperture table, and each migration aperture of migration aperture table is corresponding attenuation coefficient array in the attenuation coefficient table.The attenuation coefficient table also is a three-dimensional matrice, the corresponding offset distance of dimension, and another dimension counterparty parallactic angle, the third dimension is sequence c
j(j=1, n
i) increase to arrange with i, i=1 wherein, k, k are first element values that same offset distance and position angle are located in the migration aperture table.When picking up attenuation coefficient, offset distance and position angle (its sine value) of seismic trace rounded by sampling interval, directly in table, pick up.The position angle of seismic trace and offset distance have determined the migration aperture table, have also determined the attenuation coefficient table, and both deposit and pick up with quadrat method at employing.
The weight coefficient that described employing obtains based on the deconvolution image-forming condition of one way wave equation and depth shift squints, and crest meter realizes at last like this: definition
(x in the formula
s, y
s) be the shot point coordinate of seismic trace, (x
r, y
r) be the geophone station coordinate, (x, y T) are the level and time depth (during outward journey) coordinate of imaging point, V
RmsBe the migration velocity at imaging point place, then shot point to imaging point arrives the seimic travel time of geophone station again
Weight coefficient does
Seismic trace is used FFT (FFT), in effective band, each frequency component is taken advantage of-j ω, wherein j is a unit imaginary number, and ω is a frequency; Making the time-sampling of seismic trace record seismic signal is Δ t, in frequency field upper frequency limit is taken as 2/ Δ t, will be changed to zero above the value of actual frequency upper limit segment; Use FFT and do inverse fourier transform, the time-sampling of the seismic signal that obtain this moment becomes 0.25 Δ t; 4 τ/Δ t-0.001 rounded obtain integer i; Value by 2 of i on the time domain data after the conversion and i+1 is done linear interpolation; Obtain the value a at time τ place, the amplitude that then squints is taken as
Describe and to know through such scheme; Core of the present invention has 3 points: the one, analyze and the time migration of research three-dimensional prestack from the three-dimensional prestack depth migration method; Based on the one way ripple theoretical with surely put principle mutually; Obtain seimic travel time and amplitude in the three-dimensional non-uniform dielectric,, obtained protecting the weight coefficient of the width of cloth then based on the deconvolution image-forming condition of pre-stack depth migration; The 2nd, from the inclination angle analysis of 3-D migration impulse response, provide by along survey line and vertical line direction two maximum imaging angles (structure inclination angle) decision the time become the 3-D migration aperture; The 3rd, through the migration aperture edge being introduced banded decay, compacting skew noise automatically.Its concrete realization principle is following:
1. seimic travel time, amplitude and weight coefficient
Suppose that three-dimensional nonhomogeneous media can be approximately layered medium, in wave number one frequency field, the degree of depth continuation of the seismic wave field of single shot point or geophone station can use whilst on tour (time depth) to be expressed as:
Δ T in the formula
iBe the thickness that each layer medium expressed during with outward journey, n is the number of plies that a certain degree of depth comprises,
Be the time depth (during outward journey) of this actual grade, v
iBe the speed of each layer medium, p
xAnd p
yBe ray parameter, j is a unit imaginary number, and ω is a frequency, and f (ω) is the Fourier transform of the time-domain signal of shot point or geophone station.
(1) but in the formula in the exponential term of exponential function add up the part approximate expression be:
Taylor is on (2) formula both sides launches, and be truncated to second, can get:
V
RmsPromptly be the needed migration velocity of pre-stack time migration, i.e. root-mean-square velocity.
With (2) formula substitution (1) formula and do the space inverse fourier transform, the wave field of space one frequency field is:
Formula (4) is available surely put mutually principle try to achieve progressive separate into:
In the formula (5), definition
And
Yes
and
zero.Try to achieve
and
substitution (5) formula, in the time of can solving the walking of seismic event in the three-dimensional nonhomogeneous media and amplitude do
In the formula
Since based on phase-shift method when surely putting principle mutually and obtained walking more accurately and amplitude, protect width of cloth imaging with regard to the deconvolution image-forming condition of Wave equation depth migration capable of using.If establish focus is pulse between a period of time, then to a seismic trace, imaging results is arranged:
F ' is the first order derivative of the corresponding time-domain function of f (ω) (t) in the formula, A
sAnd t
sBe shot point (x
s, y
s, 0) and (amplitude T) can be through calculating when walking for x, y to imaging point
Try to achieve the V in the formula by (6) and (7) formula
RmsIt is the migration velocity of imaging point; A
rAnd t
rBe by the amplitude of the geophone station of trying to achieve to imaging point with quadrat method when walking.When not needing to calculate away respectively in the practical application and amplitude, directly calculate
And then calculate when always walking and weight coefficient does
Formula (8) image-forming condition shows: to arbitrary seismic trace and arbitrary imaging point, seismic signal is asked first order derivative, according to the migration velocity at imaging point place, calculate t when always walking
r+ t
sWith weight coefficient A
r/ A
s, on the first order derivative of seismic signal, pick up t
s+ t
rValue constantly also is multiplied by weight coefficient, promptly obtains the skew amplitude of this seismic trace at this imaging point, and adding up of all seismic trace skew amplitudes promptly obtains the result of migration imaging.
Be different from the current methods of ignoring weight coefficient, the weight coefficient of (8) formula has guaranteed correctly to have compensated in the migration imaging geometry diffusional effect of seismic wave propagation.
(8) F ' in the formula is a discrete function (t), picks up its t
s+ t
rValue constantly generally needs the applying interpolation technology.For the calculated amount that reduces interpolation and improve interpolation precision, the present invention adopts the resampling interpolation technique: in frequency field with maximum frequency f
Max=0.5/ Δ t enlarges 4 times, is changed to zero to the frequency component of the frequency part that increases, and does inverse fourier transform with FFT then; This moment time-sampling become 0.25 Δ t, under this sampling to 4 (t
s+ t
r)/Δ t rounds, and does linear interpolation by adjacent two point values, can obtain F ' (t very accurately
s+
Tr) value.
When the time depth T that (8) uses in the formula was outward journey, for consistent with the existing imaging results of expressing with TWT, the imaging results of (8) formula can be used I, and (x, y 2T) expressed.
2. the three-D migration aperture that becomes the time
As far as single seismic trace calculations of offset, the calculating of (8) formula is at great imaging space, and promptly (x, y T) carry out in the scope, are exactly the migration aperture of calculations of offset.Migration aperture is important to migration before stack: the calculations of offset amount can be reduced in less aperture, can not be but exist to the risk of the correct imaging of steep dip structure, and excessive aperture has been brought skew noise and bigger calculated amount again.For this reason, we have invented and have utilized the inclination maximum that is configured in both direction, the method in the three-D migration aperture that becomes when self-adaptation is confirmed.
The time three-D migration aperture that becomes be that (x y) goes up initial imaging time T with each horizontal coordinate point
bT is described
bGreater than the imaging maximum time the degree of depth the zone be exactly the zone that need not form images.
Migration aperture of the present invention is relevant with the offset distance and the position angle of seismic trace.At first by two inclination maximums of structure along survey line and vertical line direction, θ
xAnd θ
y, make up following inequation group:
In the formula
Be the position angle of this seismic trace, satisfy the minimum of alpha and the β of (9) and (10) formula, guarantee that exactly two inclination maximums are θ
xAnd θ
yStructure when being able to correctly form images, along two maximum imaging angles of azimuth direction and vertical orientations angular direction.
To see the three-D migration impulse response of single seismic trace as coordinate axis along azimuth direction and vertical orientations angular direction, can regard three-D migration impulse response as one group of ellipse is one group of ellipsoid that the axle rotation forms with the azimuth direction.Be considered to the migration aperture at picture (i.e. structure) inclination angle, be exactly only the inclination angle less than the ellipsoid part of inclination maximum zoning as (8) formula.To different time degree of depth T; Excised three-D migration impulse response greater than the inclination maximum part; Can regard as with the azimuth direction is that long axis direction, vertical orientations angular direction are short-axis direction, the elliptical area of center on the central point
of this seismic trace; Utilize following geometric relationship: major axis is that shift pulse response ellipsoid is at a string along ellipse on the azimuth direction tangent plane; This string and oval intersection point place, the inclination angle of oval tangent line is α; Can try to achieve the length of string.In like manner can try to achieve the length of the string of circle on the tangent plane of vertical orientations angular direction, the inclination angle of the tangent line that this moment is round is β.The long and short semiaxis that can try to achieve the zoning elliptical area like this is respectively
is the average root-mean-square speed on the different time degree of depth in the formula; It is half offset distance of seismic trace that parameter lambda is defined as
, when time depth T is outward journey.
To time depth T and T+ Δ T (Δ T is the sampling interval that forms images on the time depth direction); Calculate corresponding long and short semiaxis respectively; These two oval bands that form be exactly initial imaging time be the imaging region of T; To the time depth sampled point circulation (greater than this time depth, migration aperture will not change in time) less than the fixed time degree of depth, the three-D migration aperture that becomes in the time of just can obtaining.
For consistent with the imaging results of expressing with TWT, the initial imaging time in the migration aperture is expressed with 2Tb.
3. skew noise compacting
Because limiting factors such as long, the effective migration aperture of limited collection cable; The accumulation calculating of skew will produce the skew noise; The present invention adopts the input channel imaging mode, introduces the decay along the migration aperture edge before adding up through the skew amplitude at single seismic trace, the noise of realizing squinting compacting.The imaging space of single seismic trace is corresponding with migration aperture, and the normal direction that attenuation coefficient can be respectively increases the oval border of the direction and the different time degree of depth along T applies, and concrete grammar is following:
1) the definition three-dimensional matrice (x, y, T), to each horizontal coordinate point of imaging region, according to (the x that obtains from migration aperture
i, y
i) point initial imaging time T
b, design factor
And deposit it in correspondence (x
i, y
i, T) in the position, n is counting of comprising of the marginarium of artificial definition in the formula, generally is taken as 20, Δ T is the sampling of time depth direction imaging;
2) to each time depth T; Calculate a and b by (11) and (12) formula; By the horizontal direction imaging sampling interval Δ y of maximum, calculate
and then calculating again
where
is the azimuth seismic trace.(x to f>1
i, y
i, T
i), if corresponding (x
i, y
i, T) stored c in the position
1, do not do operation, otherwise calculate
With c
2Deposit in the corresponding position;
The coefficient of 3) step 2) trying to achieve is a space distribution, if the coefficient that will not have assignment point is as 1.0, it is level and smooth to be space Gauss to these coefficients, only keeps the element less than 1.0 after smoothly, forms attenuation coefficient by these elements.
To single seismic trace, to the I (x that calculates by (8) formula in the migration aperture
i, y
i, 2T),, take advantage of with attenuation coefficient that (x, y 2T), have promptly realized the edge decay, and all so single seismic trace migration result add up mutually, can suppress the skew noise effectively in skew amplitude I in the migration aperture marginarium.
Show through technique scheme; The invention has the beneficial effects as follows: 1, the inventive method can be applicable to the 3-D seismics data of many towing cables or many survey line records; Through trying to achieve the weight coefficient of protecting the width of cloth; Can in imaging process, correctly compensate the geometry diffusional effect of seismic wave propagation, obtain protecting the CRP gather of the width of cloth; Can be according to two inclination maximums of structure along survey line and vertical line direction, the three-D migration aperture that becomes when self-adaptation is confirmed, thus can reduce skew noise and the counting yield that improves offset method; Can in imaging process, suppress the skew noise automatically, obtain migration imaging result than high s/n ratio.2, the inventive method is suitable for the Parallel Implementation of cluster computer.The migrated image that the inventive method obtains can be indicated underground structure form, structure fracture location and stratum depositional pattern, to confirming favourable life, oil-bearing structure and favourable oil and gas reservoir significant application value is arranged.Guarantor's width of cloth common reflection point (CRP) road collection that this method generates can be served oil gas and fluid detection technology such as prestack inversion better.
Embodiment
1: three repair and maintenance width of cloth of embodiment prestack time migration method is an example to zone, Bohai Sea JZ32-4 block, is specially following steps:
(1) writes down the seismic reflection signals of artificial epicenter excitation with many towing cables or many surveys line, record on the tape.Specifically be, obtain the seismic reflection signals that man-made explosion excites, record on the tape with the marine streamer acquisition mode; Every big gun has 3 towing cables, between towing cable apart from 100m, the spacing of wave detector on the towing cable (being track pitch) 3.125m, the record duration 5s of seismic signal, time-sampling 1ms, every big gun contains 3 * 1152=3456 road altogether; Excite and write down 2745 big guns altogether.Fig. 1 is the typical single shot record after the conventional noise compression process.
(2) read seismic data from tape, the pre-stack seismic data is carried out conventional compacting noise treatment.To several common midpoints position, extract CMP gather, to conventional NMO correction (NMO) velocity pick of these road collection dos, the gained result is done laterally on average, as initial, laterally uniform migration velocity field.The pre-stack seismic data is pressed offset distance ordering, is divided into groups, with not on the same group seismic data be placed on the various computing node of cluster computer.Specific as follows; To be designated as
based on the laterally uniform migration velocity field that the NMO velocity pick obtains; The pre-stack seismic data 70 groups have been divided into; Every group of seismic data only comprises the seismic trace in the limited offset distance scope, and the 1st group of offset distance variation range is 150-250m, and the 2nd group is 250-350m; By that analogy, the 70th group is 3500-3600m; The data volume of respectively organizing seismic data that grouping will make as far as possible is near identical.
(3) according to the inclination maximum of the offset distance of seismic data, azimuthal variation range and underground structure and initial, horizontal migration velocity field (i.e.
) uniformly, the migration aperture table in the three-D migration aperture that becomes during structure.Specific as follows, Calculation of Three Dimensional migration aperture season is along the inclination maximum θ of side line directional structure vectorical structure
x=50 degree, the inclination maximum θ of vertical side line directional structure vectorical structure
y=30 degree, the horizontal coordinate sampling of migration imaging is Δ x=6.25m and Δ y=25m, time depth (during outward journey) sampling Δ T=0.5ms; The migration aperture table is pressed offset distance from 150 to 3600m, and position angle sine value from 0 to 0.8 calculates, and the sampling interval of offset distance and position angle sine value is respectively 100m and 0.14 in the table.When picking up migration aperture, the offset distance and the position angle sine value of seismic trace rounded by 100m and 0.14, directly in table, pick up.Fig. 2 is that offset distance is that 200m, position angle sine value are that 0.14 o'clock migration aperture is at initial imaging time (TWT) curve on survey line and vertical survey line tangent plane.
With one group of data (x
i, y
i, 2T
i) (i=1 k) describes the three-D migration aperture, and wherein k is total discrete counting in effective imaging region (imaging space is in the projection of surface level), (x
i, y
i) be that two horizontal coordinates of certain discrete point and seismic trace central point are poor in the imaging region, Ti is then to be the initial imaging time (during outward journey) of this calculations of offset; Making the shot point coordinate of seismic trace is (x
s, y
s), the geophone station coordinate is (x
r, y
r), the migration aperture that becomes during introducing is exactly in the calculations of offset of this seismic trace, is (x to horizontal coordinate in the imaging region
i+ 0.5 (x
s+ x
r), y
i+ 0.5 (y
s+ y
r)) point (i=1, k), from T=T
iBegin to carry out calculations of offset.
By two inclination maximums of structure along survey line and vertical line direction, θ
xAnd θ
y, make up following inequation group:
In the formula
Be the position angle of seismic trace, satisfy the minimum angles α and the β of last two formulas, guarantee that exactly inclination maximum is θ
xAnd θ
yStructure when being able to correctly form images, along two maximum imaging angles of azimuth direction and vertical orientations angular direction.
To different time degree of depth T; It is that long axis direction, vertical orientations angular direction are short-axis direction, the elliptical area of center on the central point of seismic trace that imaging region can be regarded as with the azimuth direction, by α, β and initial, laterally migration velocity field
can be tried to achieve its long and short semiaxis and is respectively uniformly
To be defined as
be half offset distance of seismic trace to parameter lambda in the formula, when time depth T is outward journey.
To time depth T and T+ Δ T (Δ T is the sampling interval that forms images on the time depth direction); Calculate corresponding long and short semiaxis respectively; Two ellipses that make this two group leader, minor semi-axis form have common central point (0; 0), these two oval bands that form be exactly initial imaging time be the imaging region of T; With the discrete point coordinate x=m Δ x that drops in this band, it is right that y=n Δ y and T constitute data, and wherein Δ x and Δ y are the horizontal coordinate samplings of migration imaging.To time depth sampled point circulation (greater than this time depth, migration aperture will not change in time), to the whole (x that obtain less than the fixed time degree of depth; Y; The size of 2T) pressing x and y respectively sorts, and adds up the k that always counts, the three-D migration aperture that becomes in the time of just can obtaining.
The position angle and the offset distance of seismic trace have determined migration aperture; Need position angle and offset distance variation range according to whole seismic traces; (be the minimizing memory space according to angle and offset distance sampling interval; Desirable bigger spacing), calculates the three-D migration aperture that changes with position angle and offset distance, make up the migration aperture table; The migration aperture table is a three-dimensional matrice, the corresponding offset distance of dimension, and another dimension counterparty parallactic angle, the third dimension are three coordinate data (x in three-D migration aperture
i, y
i, 2T
i) and n
iTotally 4 numbers increase arrangement, n here with i
iIt is one group of integer of corresponding attenuation coefficient; First element of this one dimension is the k that always counts of this migration aperture.That in fact, deposit in this matrix is coordinate (x
i, y
i, 2T
i) corresponding discrete point sequence number, so the migration aperture table is an INTEGER MATRICES.
Azimuthal sampling interval is not to the actual angle value, but to azimuthal sine value.When picking up migration aperture, examine sine value, offset distance and position angle sine value are rounded by sampling interval, directly in table, pick up in the coordinate difference computer azimuth angle of vertical side line direction by the offset distance and the big gun of seismic trace.
(4) confirm corresponding compacting skew noise alleviation coefficient table according to the migration aperture table.Specifically be, get the n=40 that counts that the marginarium comprises, time depth (during outward journey) sampling Δ T=0.5ms, maximum horizontal direction imaging sampling interval Δ y=25m; The attenuation coefficient table is pressed offset distance from 150 to 3600m, and position angle sine value from 0 to 0.8 calculates, and the sampling interval of offset distance and position angle sine value is respectively 100m and 0.14 in the table.When picking up attenuation coefficient, the offset distance and the position angle sine value of seismic trace rounded by 100m and 0.14, directly in table, pick up.
The attenuation coefficient array that each migration aperture is corresponding is obtained by following steps:
1) the definition three-dimensional matrice (x, y, T), according to (the x that obtains from migration aperture
i, y
i) point initial imaging time T
b, design factor
And deposit it in correspondence (x
i, y
i, T) in the position, n is counting of comprising of the marginarium of artificial definition in the formula, Δ T is the sampling interval (during outward journey) of time depth direction imaging, to all k dot cycles in the migration aperture;
2) to each time depth T, by computes a and b:
Each parameter identical when calculating migration aperture in the formula.By the horizontal direction imaging sampling interval Δ y of maximum, calculating
and then calculate
wherein,
is the position angle of seismic trace.(x to f>1
i, y
i, T), if corresponding (x
i, y
i, T) stored c in the position
1, do not do operation, otherwise calculate
With c
2Deposit in the corresponding position;
3) (each element T) circulates, and will in preceding two steps, not have the some assignment 1.0 of assignment for x, y to three-dimensional matrice; It is level and smooth that matrix element is space Gauss, only keeps the element less than 1.0 after level and smooth, forms the attenuation coefficient array by these elements.Every group of (x of migration aperture array
i, y
i, 2T
i) corresponding sequence c in the attenuation coefficient array
j(j=1, n
i), n wherein
iBe matrix (x, y, T) in (x
i, y
i) locate, along T
bPlay element number, c less than 1.0
jBe these less than 1.0 element value, as discussing in the relevant migration aperture table, n
iTo leave in the migration aperture table.
The attenuation coefficient table is corresponding with the migration aperture table, and each migration aperture of migration aperture table is corresponding attenuation coefficient array in the attenuation coefficient table.The attenuation coefficient table also is a three-dimensional matrice, the corresponding offset distance of dimension, and another dimension counterparty parallactic angle, the third dimension is sequence c
j(j=1, n
i) increase to arrange with i, i=1 wherein, k, k are first element values that same offset distance and position angle are located in the migration aperture table.When picking up attenuation coefficient, offset distance and position angle (its sine value) of seismic trace rounded by sampling interval, directly in table, pick up.The position angle of seismic trace and offset distance have determined the migration aperture table, have also determined the attenuation coefficient table, and both deposit and pick up with quadrat method at employing.
(5) weight coefficient that adopts the deconvolution image-forming condition based on one way wave equation and depth shift the to obtain amplitude calculating of squinting.
Definition
(x in the formula
s, y
s) be the shot point coordinate of seismic trace, (x
r, y
r) be the geophone station coordinate, (x, y T) are the level and time depth (during outward journey) coordinate of imaging point, V
RmsBe the migration velocity at imaging point place, then shot point to imaging point arrives the seimic travel time of geophone station again
Weight coefficient does
Seismic trace is used FFT (FFT), in effective band, each frequency component is taken advantage of-j ω, wherein j is a unit imaginary number, and ω is a frequency; Make seismic trace write down between drawing of seismic signal and be sampled as Δ t, upper frequency limit is taken as 2/ Δ t, will be changed to zero above the value of actual frequency upper limit segment in frequency field; Use FFT and do inverse fourier transform, the time-sampling of the seismic signal that obtain this moment becomes 0.25 Δ t; 4 τ/Δ t-0.001 rounded obtain integer i; Value by 2 of i on the time domain data after the conversion and i+1 is done linear interpolation; Obtain the value a at time τ place, the amplitude that then squints is taken as
(6) at each computing node, to whole seismic trace circulations.To each seismic trace, confirm the initial imaging time of each discrete point calculations of offset in effective imaging region and the zone by the migration aperture table, begin to calculate the skew amplitude of each imaging point from initial imaging time.To the imaging point in the marginarium that is in the three-D migration aperture, look into the attenuation coefficient table pick up attenuation coefficient and multiply by the skew amplitude obtain final amplitude.Amplitude is added on the offset distance of depositing correspondence in the migration result array.Specifically be, define half offset distance separation delta h=12.5m, calculate m=(h
Max-h
Min)/Δ h+1, wherein h
MaxBe maximum half offset distance of all seismic traces on this node, h
MinBe smaller part offset distance, if three-dimensional imaging space size is (n
x, n
y, n
T), then this node array of depositing migration result is (n
x, n
y, n
T, m); Migration result will be added to the m of this array the 4th dimension according to half offset distance h of each seismic trace
0=(h-h
MinOn the)/Δ h+1.Seismic data can not distinguished independent calculating on the various computing node of cluster computer on the same group, and corresponding migration aperture table and attenuation coefficient table only need be deposited the corresponding offset distance of this group seismic data and the part of azimuth coverage.
(7) collect the calculations of offset result of each computing node, to each horizontal level of imaging region the calculations of offset result by the time degree of depth, offset distance ordering, form the CRP gather of each horizontal level.Specifically be, with the matrix (n that deposits migration result on each computing node
x, n
y, n
T, m
i) be integrated into a matrix (n
x, n
y, n
T, ∑ m
i), wherein the length of the 4th dimension is adding up of each computing node this matrix corresponding dimension length; To the identical part of offset distance of seismic data on the various computing node, will be superimposed these results, be stored in then on the corresponding position of matrix, the length of matrix the 4th dimension is also wanted corresponding minimizing.Each horizontal level (n
xΔ x, n
yΔ y) two-dimensional array (n that locates
T, ∑ m
i) promptly be the CRP gather of this horizontal level.
(8) whether straight according to the lineups in the CRP gather, in migration process, confirm the migration velocity field.Specifically be; Adopt laterally uniform migration velocity field
to carry out the calculations of offset first time; CRP gather to the skew generation; Do inverse dynamic correction according to
; Do NMO correction again and obtain new speed; This speed is done space smoothing handle, promptly obtain final migration velocity field.Fig. 3 is the migration velocity field of trying to achieve along the isogram on the survey line tangent plane.
(9) utilize the migration velocity field of confirming in the migration process to do calculations of offset; The whole CRP gathers that form are done the moving school of residue; The corresponding numerical value that obviously stretching and noise part occur is changed to zero (i.e. excision), with the numerical value stack of different offset distances, formation skew stacked section.Fig. 4 is through the typical CRP gather after the moving school of residue, stretching and the noise excision.
(10) be the profile image of underground reflective construct through the software for display stacked section data-switching that will squint, profile image will be indicated form, breakpoint, turn-off size and the interface of the underground structure reflection strength to seismic event.Fig. 5 constructs the profile image along the line direction tangent plane with the underground buried hill of JZ32-4 block, zone, the Bohai Sea that the inventive method obtains, and buried hill flank, top, tomography, stratum depositional pattern are well portrayed.