CN101957455B - Method of three-dimensional preserved-amplitude pre-stack time migration - Google Patents

Abstract

A method of three-dimensional preserved-amplitude pre-stack time migration which is applied to reflected seismic data treatment in seismic prospecting is a pre-stack migration imaging method aiming at the three-dimensional collected seismic data. The technical scheme comprises three aspects: 1) performing the migration amplitude calculation by using a weight coefficient acquired according to a one-way wave equation and the reverse convolution imaging condition of the depth migration, 2) confirming the time varying three-dimensional migration aperture according to the biggest incidence formed underground, and 3) automatically suppressing the migration voice by performing strip attenuation method on the edge of the migration aperture. The method of the invention realizes high-efficient, high SNR (Signal to Noise Ratio) and preserved-amplitude migration imaging by improving the migration aperture, providing a weight coefficient of the preserved-amplitude and suppressing the migration voice in the imaging process. The preserved-amplitude reflective spot collection generated by using the method of the invention can better serve for the oil gas and fluid detection technique for pre-stack inversion. The method of the invention has high application value to oil gas and mine resource exploration.

Description

Three repair and maintenance width of cloth prestack time migration methods

Technical field

The invention belongs to reflected seismic information processing technology field in the oil gas field seismic prospecting, relate in particular to the migration before stack imaging technique category in the seismic data processing procedure, is three repair and maintenance width of cloth prestack time migration methods of 3-D seismics data migration imaging.

Background technology

In the reflected seismic information processing procedure, the migration before stack imaging is the key link in the seismic prospecting, and prestack time migration method is a kind of important method in the migration before stack imaging.Prestack time migration method can be to one type of tomography complicacy but the speed horizontal change is not very violent tectonic structure better forms images comparatively.Compare with prestack depth migration method, except that having higher counting yield, its main advantage is only to need to use stack (root mean square) speed; Can obtain appropriate migration velocity model through modes such as velocity sweepings simply like this, avoid a main difficulty using prestack depth migration method to face: the speed modeling.Therefore, prestack time migration method has become the gordian technique of field of seismic exploration widespread use.

The factor that influences the pre-stack time migration imaging effect comprises: the migration aperture of migration velocity, seimic travel time computing method, calculations of offset, weight coefficient, migration algorithm realization flow when calculating the skew amplitude.Can calculate with migration velocity when walking have determined reflection wave correctly playback jointly; Migration aperture and application mode thereof have determined the calculated amount of skew noise and migration algorithm; Weight coefficient has determined the physical parameter that can the imaging amplitude the underground actual interface of correct response to change, and the migration algorithm realization flow has material impact to counting yield and the storage demand of skew.As far as offset method, imaging effect, counting yield and storage demand are the important indicators of estimating offset method.

Migration aperture is important to pre-stack time migration, and migration aperture has more vital role to compacting skew noise and reduction calculations of offset amount preferably.Little migration aperture can reduce the calculations of offset amount, can not be to the risk of the correct imaging of steep dip structure but exist; Excessive aperture has been brought skew noise and bigger calculated amount again.As far as the three-dimensional prestack time migration method, because image data is having different spatial sampling rates along survey line and vertical line direction, choosing of migration aperture is complicated more than two-dimensional migration.

Existing prestack time migration method is the basis with the double square root equation, and in fact its image-forming condition has used the dependent imaging condition of pre-stack depth migration, does not compensate the geometry diffusional effect of seismic wave propagation when therefore forming images.Protect for obtaining common reflection point (CRP) the road collection of the width of cloth, utilize prestack information to carry out oil gas and fluid detection, need development to protect the width of cloth, be i.e. the prestack time migration method that can correct response subsurface interface physical parameter changes of amplitude to serve prestack inversion etc.

Three-dimensional observation has become the main observed pattern of marine and land seismic prospecting, and three-dimensional observation and three-D migration can be to the correct imaging of steep dip structure, the problems that lateral reflection is difficult to correct playback when having avoided two-dimensional observation and skew.Along with the development of collection and computing equipment, what three-dimensional acquisition and migration processing had become realizes more easily, and it is technological therefore must to develop corresponding three-dimensional prestack time migration.

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

$a=\frac{{\stackrel{\‾}{V}}_{\mathrm{rms}}\left(T\right)\·T}{\sin \left(2\mathrm{\α}\right)}\sqrt{1+{\cos }^2\left(2\mathrm{\α}\right)+{\mathrm{\λ}}^2{\sin }^2\left(2\mathrm{\α}\right)-2\cos \left(2\mathrm{\α}\right){[1+{\mathrm{\λ}}^2{\sin }^2\left(2\mathrm{\α}\right)]}^{1/2}}$

$b={\stackrel{\‾}{V}}_{\mathrm{rms}}\left(T\right)\·T\tan \mathrm{\β}$

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:

$a=\frac{{\stackrel{\‾}{V}}_{\mathrm{rms}}\left(T\right)\·T}{\sin \left(2\mathrm{\α}\right)}\sqrt{1+{\cos }^2\left(2\mathrm{\α}\right)+{\mathrm{\λ}}^2{\sin }^2\left(2\mathrm{\α}\right)-2\cos \left(2\mathrm{\α}\right){[1+{\mathrm{\λ}}^2{\sin }^2\left(2\mathrm{\α}\right)]}^{1/2}}$

$b={\stackrel{\‾}{V}}_{\mathrm{rms}}\left(T\right)\·T\tan \mathrm{\β}$

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 ₂Deposit 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

$p_s=\frac{{(x-x_s)}^2+{(y-y_s)}^2}{V_{\mathrm{rms}}^2{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="HSA00000279585500021.png,HSA00000279585500051.png"\; state="\{\{state\}\}">T</figure-callout>}}^2},$ $p_r=\frac{{(x-x_r)}^2+{(y-y_r)}^2}{V_{\mathrm{rms}}^2{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="HSA00000279585500021.png,HSA00000279585500051.png"\; state="\{\{state\}\}">T</figure-callout>}}^2}$

(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:

$\tilde{P}(p_x,p_y,\mathrm{\&omega;},T=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;T}}_i)=f\left(\mathrm{\&omega;}\right)\exp (-\mathrm{j\&omega;}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;T}}_i\sqrt{1-v_i^2(p_x^2+p_y^2)})---\left(1\right)$

Δ 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:

$\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left({\mathrm{\&Delta;T}}_i\sqrt{1-v_i^2(p_x^2+p_y^2)}\right)\&ap;\sqrt{1-V_{\mathrm{rms}}^2(p_x^2+p_y^2)}\&CenterDot;\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;T}}_i\right)---\left(2\right)$

Taylor is on (2) formula both sides launches, and be truncated to second, can get:

$V_{\mathrm{<figure-callout\; id="2"\; label="rms"\; filenames="HSA00000279585500021.png,HSA00000279585500051.png"\; state="\{\{state\}\}">rms</figure-callout>}}^2=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{v}_i^2{\mathrm{\&Delta;T}}_i}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;T}}_i}---\left(3\right)$

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:

$P(x,y,\mathrm{\&omega;},T)=\frac{\mathrm{\&omega;}}{{4\mathrm{\&pi;}}^2}\&Integral;\&Integral;f\left(\mathrm{\&omega;}\right)\exp [-\mathrm{j\&omega;}(T\sqrt{1-V_{\mathrm{rms}}^2(p_x^2+p_y^2)}-p_{x}x-p_{y}y)]{\mathrm{dp}}_x{\mathrm{dp}}_y---\left(4\right)$

Formula (4) is available surely put mutually principle try to achieve progressive separate into:

$P(x,y,\mathrm{\&omega;},T)=f\left(\mathrm{\&omega;}\right)\frac{\mathrm{\&omega;}}{2\mathrm{\&pi;}}\exp (-j\frac{\mathrm{\&pi;}}{2}){\left|Q(p_x^0,p_y^0)\right|}^{-1/2}\exp [-j\mathrm{\&omega;\&phi;}(p_x^0,p_y^0)]---\left(5\right)$

In the formula (5), definition

$\mathrm{\&phi;}(p_x,p_y)=T\sqrt{1-V_{\mathrm{rms}}^2(p_x^2+p_y^2)}-p_{x}x-p_{y}y$

$Q(p_x,p_y)=\left|\begin{array}{cc}\frac{{\&PartialD;}^2\mathrm{\&phi;}}{{\&PartialD;p}_x^2}& \frac{{\&PartialD;}^2\mathrm{\&phi;}}{{\&PartialD;p}_x{\&PartialD;p}_y}\\ \frac{{\&PartialD;}^2\mathrm{\&phi;}}{{\&PartialD;p}_x{\&PartialD;p}_y}& \frac{{\&PartialD;}^2\mathrm{\&phi;}}{{\&PartialD;p}_y^2}\end{array}\right|$

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

$t=\mathrm{<figure-callout\; id="1"\; label="T"\; filenames="HSA00000279585500021.png"\; state="\{\{state\}\}">T</figure-callout>}\sqrt{1+p_0}---\left(6\right)$

$A=\frac{1}{{\mathrm{TV}}_{\mathrm{rms}}^2(1+p_0)}---\left(7\right)$

In the formula

${\mathrm{<figure-callout\; id="0"\; label="p"\; filenames="HSA00000279585500011.png,HSA00000279585500021.png"\; state="\{\{state\}\}">p</figure-callout>}}_0=\frac{x^2+y^2}{V_{\mathrm{rms}}^2{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="HSA00000279585500021.png,HSA00000279585500051.png"\; state="\{\{state\}\}">T</figure-callout>}}^2},$

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:

$I(x,y,T)=\frac{A_r}{A_s}\&Integral;f\left(\mathrm{\&omega;}\right)\mathrm{\&omega;exp}(-j\frac{\mathrm{\&pi;}}{2})\exp (-\mathrm{j\&omega;}(t_s+t_r))\mathrm{d\&omega;}---\left(8\right)$

$=(A_r/A_s)F^{\&prime;}(t_s+t_r)$

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

$p_s=\frac{{(x-x_s)}^2+{(y-y_s)}^2}{V_{\mathrm{rms}}^2{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="HSA00000279585500021.png,HSA00000279585500051.png"\; state="\{\{state\}\}">T</figure-callout>}}^2},$ $p_r=\frac{{(x-x_r)}^2+{(y-y_r)}^2}{V_{\mathrm{rms}}^2{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="HSA00000279585500021.png,HSA00000279585500051.png"\; state="\{\{state\}\}">T</figure-callout>}}^2}$

And then calculate when always walking and weight coefficient does

$t_s+t_r=(\sqrt{1+p_s}+\sqrt{1+p_r})T$

$\frac{A_r}{A_s}=\frac{1+p_s}{{1+p}_r}$

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

$a=\frac{{\stackrel{\&OverBar;}{V}}_{\mathrm{rms}}\left(T\right)\&CenterDot;T}{\sin \left(2\mathrm{\&alpha;}\right)}\sqrt{1+{\cos }^2\left(2\mathrm{\&alpha;}\right)+{\mathrm{\&lambda;}}^2{\sin }^2\left(2\mathrm{\&alpha;}\right)-2\cos \left(2\mathrm{\&alpha;}\right){[1+{\mathrm{\&lambda;}}^2{\sin }^2\left(2\mathrm{\&alpha;}\right)]}^{1/2}}---\left(11\right)$

$b={\stackrel{\&OverBar;}{V}}_{\mathrm{rms}}\left(T\right)\&CenterDot;T\tan \mathrm{\&beta;}---\left(12\right)$

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 ₁, do not do operation, otherwise calculate

With c ₂Deposit 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.

Description of drawings

Fig. 1 is the typical single shot record after the conventional noise compression process in JZ32-4 area, zone, the Bohai Sea, and collection is 3 towing cables of every big gun;

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;

Fig. 3 is the migration velocity field along the isogram on the survey line tangent plane, and numeral is the migration velocity value among the figure;

Fig. 4 is through the typical CRP gather after the moving school of residue, stretching and the noise excision;

Fig. 5 is the profile image of zone, Bohai Sea JZ32-4 block underground buried hill structural offset result along the line direction tangent plane.

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

$a=\frac{{\stackrel{\&OverBar;}{V}}_{\mathrm{rms}}\left(T\right)\&CenterDot;T}{\sin \left(2\mathrm{\&alpha;}\right)}\sqrt{1+{\cos }^2\left(2\mathrm{\&alpha;}\right)+{\mathrm{\&lambda;}}^2{\sin }^2\left(2\mathrm{\&alpha;}\right)-2\cos \left(2\mathrm{\&alpha;}\right){[1+{\mathrm{\&lambda;}}^2{\sin }^2\left(2\mathrm{\&alpha;}\right)]}^{1/2}}$

$b={\stackrel{\&OverBar;}{V}}_{\mathrm{rms}}\left(T\right)\&CenterDot;T\tan \mathrm{\&beta;}$

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:

$a=\frac{{\stackrel{\&OverBar;}{V}}_{\mathrm{rms}}\left(T\right)\&CenterDot;T}{\sin \left(2\mathrm{\&alpha;}\right)}\sqrt{1+{\cos }^2\left(2\mathrm{\&alpha;}\right)+{\mathrm{\&lambda;}}^2{\sin }^2\left(2\mathrm{\&alpha;}\right)-2\cos \left(2\mathrm{\&alpha;}\right){[1+{\mathrm{\&lambda;}}^2{\sin }^2\left(2\mathrm{\&alpha;}\right)]}^{1/2}}$

$b={\stackrel{\&OverBar;}{V}}_{\mathrm{rms}}\left(T\right)\&CenterDot;T\tan \mathrm{\&beta;}$

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 ₁, do not do operation, otherwise calculate With c ₂Deposit 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

$p_s=\frac{{(x-x_s)}^2+{(y-y_s)}^2}{V_{\mathrm{rms}}^2{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="HSA00000279585500021.png,HSA00000279585500051.png"\; state="\{\{state\}\}">T</figure-callout>}}^2},$ $p_r=\frac{{(x-x_r)}^2+{(y-y_r)}^2}{V_{\mathrm{rms}}^2{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="HSA00000279585500021.png,HSA00000279585500051.png"\; state="\{\{state\}\}">T</figure-callout>}}^2}$

(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 ₀=(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.

Claims (1)

1. repair and maintenance width of cloth prestack time migration method is characterized in that adopting following steps: 1) write 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 made conventional NMO velocity pick up, 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 level; 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, and the numerical value stack with different offset distances forms the stacked section that squints; 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 the described offset distance according to seismic data of step 3), azimuthal variation range and underground structure and initial, laterally uniform migration velocity field, the migration aperture table in the three-D migration aperture that becomes during structure is achieved in that with one group of data (x _i, y _i, 2T _i) the three-D migration aperture described, i=1 wherein ..., k, k are discrete counting total in effective imaging region, effectively imaging region is the projection of imaging space at 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 of this discrete point 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, from T=T _iBegin to carry out calculations of offset, i=1 wherein ..., k; By the inclination maximum θ of structure along line direction _xWith inclination maximum θ along vertical line direction _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 the maximum imaging angle [alpha] of azimuth direction with along the maximum imaging angle beta of 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 is 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

$a=\frac{{\stackrel{\&OverBar;}{V}}_{\mathrm{rms}}\left(T\right)\&CenterDot;T}{\sin \left(2\mathrm{\&alpha;}\right)}\sqrt{1+{\cos }^2\left(2\mathrm{\&alpha;}\right)+{\mathrm{\&lambda;}}^2{\sin }^2\left(2\mathrm{\&alpha;}\right)-2\cos \left(2\mathrm{\&alpha;}\right){[1+{\mathrm{\&lambda;}}^2{\sin }^2\left(2\mathrm{\&alpha;}\right)]}^{1/2}}$

$b={\stackrel{\&OverBar;}{V}}_{\mathrm{rms}}\left(T\right)\&CenterDot;T\tan \mathrm{\&beta;}$

Parameter lambda is defined as in the formula $\mathrm{\&lambda;}=h/({\stackrel{\&OverBar;}{V}}_{\mathrm{Rms}}\left(T\right)\&CenterDot;T),$ $h=0.5\sqrt{{(x_s-x_r)}^2+{(y_s-y_r)}^2}$ Be half offset distance of seismic trace, 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 less than the fixed time degree of depth; To obtain whole (x, y 2T) press the size ordering of x and y respectively; Add up the k that always counts, the three-D migration aperture that must then become;

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, according to angle and offset distance sampling interval, calculate 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 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;

Azimuthal sampling interval is to azimuthal sine value; When picking up migration aperture; Offset distance by seismic trace is examined the sine value in the coordinate difference computer azimuth angle of vertical side line direction with big gun, and offset distance and position angle sine value are rounded by sampling interval, directly in the migration aperture table, picks up;

Step 4) is described confirms that according to the migration aperture table corresponding compacting skew noise alleviation coefficient 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

T＜T wherein _b+ n Δ T, 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, the sampling interval of time depth direction imaging when Δ T is outward journey is to all k dot cycles in the migration aperture;

2) to each time depth T, by computes a and b:

$a=\frac{{\stackrel{\&OverBar;}{V}}_{\mathrm{rms}}\left(T\right)\&CenterDot;T}{\sin \left(2\mathrm{\&alpha;}\right)}\sqrt{1+{\cos }^2\left(2\mathrm{\&alpha;}\right)+{\mathrm{\&lambda;}}^2{\sin }^2\left(2\mathrm{\&alpha;}\right)-2\cos \left(2\mathrm{\&alpha;}\right){[1+{\mathrm{\&lambda;}}^2{\sin }^2\left(2\mathrm{\&alpha;}\right)]}^{1/2}}$

$b={\stackrel{\&OverBar;}{V}}_{\mathrm{rms}}\left(T\right)\&CenterDot;T\tan \mathrm{\&beta;}$

Each parameter identical when calculating migration aperture in the formula, the horizontal direction imaging sampling interval Δ y by maximum calculates

And then calculate

Wherein

Be the position angle of seismic trace, to (the x of f＞1 _i, y _i, T), if corresponding (x _i, y _i, T) stored c in the position ₁, do not do operation, otherwise calculate

With c ₂Deposit 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 wherein ..., n _i, n _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, n _iTo leave in the migration aperture table;

The attenuation coefficient table is corresponding with the migration aperture table; Each migration aperture of migration aperture table is corresponding attenuation coefficient array in the attenuation coefficient table, and the attenuation coefficient table is a three-dimensional matrice, the corresponding offset distance of dimension; Another dimension counterparty parallactic angle, the third dimension is sequence c _jIncrease to arrange with i, j=1 wherein ..., n _i, i=1 ..., k, k are first element values at same offset distance and position angle place in the migration aperture table; When picking up attenuation coefficient, the offset distance and the position angle sine value of seismic trace rounded by sampling interval, directly in the attenuation coefficient table, pick up;

The weight coefficient that the described employing of step 5) 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

$p_s=\frac{{(x-x_s)}^2+{(y-y_s)}^2}{V_{\mathrm{rms}}^2{T}^2},$ $p_r=\frac{{(x-x_r)}^2+{(y-y_r)}^2}{V_{\mathrm{rms}}^2{T}^2}$

(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 the time depth 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, 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

Images