CN102866421B - Identify the scattering wave Prestack Imaging method of little turn-off breakpoint - Google Patents

Abstract

Identifying the scattering wave Prestack Imaging method of little turn-off breakpoint, be applied to reflected seismic information process in seismic prospecting, is the seismic imaging method identifying subsurface geological structure micro-small fault.The method utilizes the breakpoint scattering wave of tomography, offsets to typical reflection ripple the independent imaging of breakpoint that impalpable turn-off is less than 1/4 wavelength.The method utilizes existing reflected seismic information, by compacting transmitted wave energy in migration before stack calculates, and utilizes the reversal of poles feature of breakpoint scattering wave, realizes the independent imaging of little breakpoint.Based on pre-stack time migration, by generating angle gathers to single big gun data skew and excise Fresnel zone corresponding to reflection line-ups, realize rejecting transmitted wave.The method can indicate minor fault, the breakpoint of underground structure subtly, thus can be familiar with migration and the shutoff of underground oil and gas better, has significant application value to the exploration and development of hydrocarbon resources.Significantly improve the recognition capability of method of seismic prospecting to underground deep crack and tomography.

Description

Identify the scattering wave Prestack Imaging method of little turn-off breakpoint

Technical field

The invention belongs to reflected seismic information processing technology field in seismic prospecting, relating to the prestack migration image technology category in process of seismic data processing, is a kind of offset imaging method that can identify the tectonic structure micro-small fault that conventional migration technique method is difficult to clearly distinguish.

Background technology

In oil-gas exploration, to the accurate identification of minor fault, little fracture, be understanding oil gas translo-cation system, and then identify the important step of Favorable Reservoir.Therefore, in seismic prospecting in reflected seismic information processing procedure, improve the resolution of migration imaging is an important effort target always.

Existing seismic data offset imaging method is all utilize the correct playback of subsurface interface reflection wave to be embodied as picture.Interface echo is in fact the stack result of multiple scattering point scattering wave, it is characterized in that an incident wave, and the reflection wave of its correspondence is only along a direction injection.

For improving the resolution of seismic migration imaging, develop many methods, comprised the spectral whitening deconvolution of imaging stacked section, unstable state deconvolution, Corpus--based Method hypothesis or all kinds of of well-log information and widen band technology, inverse Q filtering, viscoelasticity pre-stack depth migration etc.The core of each class methods is all the upper limit by improving seismic event frequency band, and what the reflection line-ups playbacked was become is thinner, thus improves the resolution to thin layer, tomography, gap.Although this type of effort has obtained certain success, the lifting of seismic event frequency band front end is always limited, still can not meet oil-gas exploration and the needs of exploitation to minor fault identification.

In fact, in the reflected seismic information of ground table record, except the interface echo that common seismic offset method utilizes, also contains the scattering wave that a large amount of small scale media variations, little breakpoint etc. produce.Due to these scattering waves amplitude comparatively reflection wave differed from one to two magnitudes, in common seismic migration processing, the playback result of these scattering waves is fallen into oblivion; And inappropriate process of seismic data processing has also suppressed scattering wave signal.

Target of the present invention is exactly the scattering wave utilizing little breakpoint etc. to produce, and to the independent imaging of these little breakpoints, thus identifies that conventional migration technique method is difficult to the little turn-off breakpoint distinguished.For realizing this target, interface echo imaging must be suppressed in calculations of offset.

The method of suppressing interface echo imaging rejects an interface echo in seismic data, this is some feature utilizing reflection wave in Prestack seismic data, as the Hyperbolic Feature in common midpoint gather realizes.But what in complex geological structure situation, the feature of interface echo became is more complicated, particularly when seismic data signal to noise ratio (S/N ratio) is lower, it is almost the work be difficult to that accurate perception reflex ripple lineups also correctly reject it; In addition, the transform methods such as radon transform are adopted also will to destroy originally very weak scattering wave.Therefore the present invention adopts and concentrates depressor reflex ripple in migration imaging road, realizes scattering wave imaging.

Because scattering wave amplitude and all kinds of noise are in same magnitude, even compared with some relevant noises also low magnitude.For avoiding noise to the interference of scattering wave imaging results, the present invention utilizes the scattering wave of little breakpoint to there is the feature of reversal of poles, and namely the polarity of scattering wave lineups is in the anti-raw positive and negative conversion of certain point, strengthens the imaging to little breakpoint.The present invention is positioned to carry out imaging to the scattering wave of little turn-off breakpoint, and these small breakpoints, be vital to identification of hydrocarbon translo-cation system.

Pre-stack time migration is a kind of important method in seismic migration imaging, and compared with prestack depth migration method, except having higher counting yield, its main advantage only need use superposition (root mean square) speed; Appropriate rate pattern can be obtained like this simply by modes such as velocity sweepings.Because the actual Geological Problems in large quantities run into is all that reflective construct is complicated, but the horizontal change of medium velocity is not very violent, and therefore the present invention develops scattering wave prestack migration image method based on prestack time migration method.This has avoided the main difficulty of a migration imaging: velocity modeling.This point is even more important to scattering wave imaging, because the slight error of scattering wave to speed is more responsive, need use the more accurate rate pattern of more conventional pre-stack depth migration.

Summary of the invention

The object of the invention is: provide a kind of scattering wave Prestack Imaging method identifying little turn-off breakpoint, it can suppress transmitted wave energy in the migration process of seismic data, obtains the independent imaging to little turn-off breakpoint; This method can improve the recognition capability of reflected seismic information migration imaging to micro-small fault, improves the understanding to Deep Oil And Gas Exploration translo-cation system in oil-gas exploration and development.

The technical solution used in the present invention is: the scattering wave Prestack Imaging method identifying little turn-off breakpoint, and concrete steps comprise:

(1) with the seismic reflection signals of the artificial epicenter excitation of acquisition mode record of single towing cable or survey line, be recorded on tape;

(2) read seismic signal from tape, form Prestack seismic data, conventional sound attenuation process is done to Prestack seismic data; Build a two dimensional cross-section by along line direction and depth direction, two dimensional processing methods is used to seismic data; For part common midpoint, extract common midpoint gather, conventional NMO velocity pickup is done to the common midpoint gather extracted, the reciprocal value of obtained speed is done laterally average, using the laterally uniform velocity field that obtains as initial offset speed;

(3) according to initial offset speed, carry out calculations of offset by conventional prestack time method, form CRP gather;

(4) to CRP gather, utilize initial offset speed to do inverse dynamic correction, then do normal moveout correction and obtain the speed after upgrading, do space smoothing to the reciprocal value of speed after the renewal of each CDP point, namely the level and smooth velocity field obtained is migration velocity field;

(5) according to migration velocity field, conventional pre-stack time migration is done to Prestack seismic data and calculates, the dynamic school of high-order residue and the excision that stretches are done to the CRP gather obtained, the start offset distance of record residual normal moveout and the excision that stretches; Along offset distance superposition CRP gather, obtain migration stack section;

(6) according to migration stack section and migration velocity field, determining the angle of lineups and surface level on migrated section, take counter clockwise direction as positive-angle;

(7) according to migration velocity field, big gun territory pre-stack time migration is done to Prestack seismic data, form the two-dimensional migration road collection with shot point and angle change at each CDP point;

(8) according to the angle of migration velocity field, shot point coordinate, imaging point coordinate, imaging point place skew lineups and surface level, set up the three-dimensional table determining Fresnel zone, three-dimensional table determines the angle domain Fresnel zone of each imaging point accordingly;

(9) balancing energy is done to the two-dimensional migration road collection of each CDP point, first determine different time degree of depth upper angle territory Fresnel zone, then on angular deflection road collection, excise Fresnel zone part, the remainder on Fresnel zone both sides is superposed respectively, obtain two imaging roads at each CDP point;

(10) on the basis at the lineups inclination angle that step 6 obtains, whether the scan method increased and decreased by number percent, completely eliminated reflection line-ups according to scattering wave migration stack section, determine to offset lineups inclination angle accurately;

(11) the accurate skew lineups inclination angle utilizing step 10 to obtain, repeat step 9, obtain two migration stack sections be made up of imaging road on the right side of imaging road on the left of Fresnel zone and Fresnel zone respectively, this two section applied energies equilibrium is subtracted each other, obtains only to the scattering wave migration stack section of the independent imaging of little turn-off breakpoint;

(12) by software for display, scattering wave migration stack section numerical value is converted to the profile image of the little breakpoint distribution of subsurface geological structure, this profile image is the tangent plane picture of structure along line direction, this image can indicate the impalpable turn-off of typical reflection ripple migrated image to be less than minor fault and the breakpoint of 1/4 wavelength, for being familiar with migration and the shutoff situation of underground oil and gas, determine effective oil and gas reservoir.

Described foundation migration velocity field, big gun territory pre-stack time migration is done to Prestack seismic data, the two-dimensional migration road collection formed with shot point and angle change at each CDP point is achieved in that each CDP point, define a three-dimensional array and deposit migration result, the corresponding time depth of one dimension, another dimension is shot point numbering, and the third dimension is angle, and angle is from bearing 60 degree to positive 60 degree of changes; Prestack seismic data is pressed common-shot-gather sequence; To the common-shot-gather of monolateral record, utilize big gun, the cautious principle such as mutual of exchanging, the common-shot-gather of the bilateral record of constructing virtual; By shot point coordinate, common-shot-gather is divided into groups, difference group common-shot-gather is placed on the different computing nodes of cluster computer; At each computing node, whole seismic trace is circulated; To each seismic trace, to the whole CDP dot cycles in migration aperture; To each CDP point, determine initial imaging point by migration aperture, the imaging point below this imaging point is calculated successively; First by shot point, geophone station and imaging point coordinate and migration velocity field, when determining to walk and imaging weight coefficient, and then determined to offset amplitude A by the seismologic record amplitude of this seismic trace; Then according to the offset distance of seismic trace, the horizontal coordinate of imaging point and time depth, the residual normal moveout that pickup step 5 records and the start offset distance stretching excision; The time depth of foundation offset distance and imaging point, judge whether skew amplitude A should be cut, if do not excise, then the angle calculating shot point incident wave and reflecting interface normal is

${\cos }^2\mathrm{\θ}=\frac{0.5}{\sqrt{1+p_s}\sqrt{1+p_r}}(1+0.5(p_s+p_r-\frac{H^2}{T^2{V}_{\mathrm{rms}}^2}))+0.5$

In above formula, θ is the angle of shot point incident wave and reflecting interface normal, unit degree, and T is the time depth of imaging point, namely during outward journey, and unit second, V _rmsfor imaging point (x ₀, T) and the migration velocity at place, unit meter per second, H is the offset distance of seismic trace, unit rice, p _sand p _rfor

$p_s=\frac{{(x_0-x_s)}^2}{V_{\mathrm{rms}}^2{T}^2},$ $p_r=\frac{{(x_0-x_g)}^2}{V_{\mathrm{rms}}^2{T}^2}$

Wherein x ₀, x _gand x _sthe horizontal coordinate of imaging point, geophone station and shot point respectively, unit rice;

If Δ θ is the angular interval of predefined, unit degree, according to cos ²(k Δ θ-Δ θ 2)≤cos ²θ <cos ²(k Δ θ-Δ θ/2) determine positive integer k; If the big gun belonging to this seismic trace is numbered n, if the residual normal moveout of pickup is τ, its unit is second, works as x _g>=x _s, then skew amplitude A is added on (T+ τ, n, the k Δ θ) element of the three-dimensional array depositing migration result, namely in positive-angle; Work as x _g<x _s, then skew amplitude A is added on (T+ τ, the n ,-k Δ θ) element of the three-dimensional array depositing migration result; Complete to whole seismic trace circulation, just obtain the two-dimensional migration road collection with shot point and angle change at each CDP point, this angle is the angle of shot point incident wave and reflecting interface normal.

The described angle according to migration velocity field, shot point coordinate, imaging point coordinate, imaging point place skew lineups and surface level, set up the three-dimensional table determining Fresnel zone, three-dimensional table is determined that angle domain Fresnel zone that each imaging point is corresponding is achieved in that and is made ρ be the ratio of imaging point place actual layer speed and migration velocity accordingly, α is the inclination angle of imaging point place reflecting interface, namely offseting the angle of lineups and surface level, take counter clockwise direction as positive-angle; Be that three dimensions set up three-dimensional table with sin φ, ρ and sin α, deposit two dimensionless factors in table, its value is calculated by following two formulas:

$a_1=\frac{\mathrm{\&rho;}(p_1\sin \mathrm{\&phi;}\cos \mathrm{\&beta;}+p_2\sin \mathrm{\&beta;}\cos \mathrm{\&phi;})-\sin \mathrm{\&alpha;}(\cos \mathrm{\&beta;}+\cos \mathrm{\&phi;})}{\mathrm{\&rho;}\cos \mathrm{\&phi;}\cos \mathrm{\&beta;}(p_1^2\cos \mathrm{\&phi;}-p_2^2\cos \mathrm{\&beta;})-2p_1\sin \mathrm{\&alpha;}\sin \mathrm{\&phi;}\cos \mathrm{\&beta;}+\left[2p_3\cos \mathrm{\&phi;}(\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}-\mathrm{\&rho;}{p}_2)\right]/{\cos }^2\mathrm{\&beta;}}$

$a_2=\frac{\mathrm{\&rho;}\cos \mathrm{\&phi;}\cos \mathrm{\&beta;}(\cos \mathrm{\&phi;}+\cos \mathrm{\&beta;})}{\mathrm{\&rho;}\cos \mathrm{\&phi;}\cos \mathrm{\&beta;}(p_1^2\cos \mathrm{\&phi;}-p_2^2\cos \mathrm{\&beta;})-2p_1\sin \mathrm{\&alpha;}\sin \mathrm{\&phi;}\cos \mathrm{\&beta;}+\left[2p_3\cos \mathrm{\&phi;}(\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}-\mathrm{\&rho;}{p}_2)\right]/{\cos }^2\mathrm{\&beta;}}$

In upper two formulas

$\sin \mathrm{\&beta;}=\sin \left(2\mathrm{\&alpha;}\right)\sqrt{1/p^2-{\sin }^2\mathrm{\&phi;}}-\cos \left(2\mathrm{\&alpha;}\right)\sin \mathrm{\&phi;},$ $\cos \mathrm{\&beta;}=\sqrt{1-{\sin }^2\mathrm{\&beta;}},$ $\cos \mathrm{\&phi;}=\sqrt{1-{\sin }^2\mathrm{\&phi;}}$

p ₁＝cosαcosφ+sinαsinφ/ρ，p ₂=cosαcosβ+sinαsinβ/ρ

$p_3=p_1\cos \mathrm{\&phi;}[\frac{\mathrm{\&rho;}\sin \left(2\mathrm{\&alpha;}\right)\sin \mathrm{\&phi;}\cos \mathrm{\&phi;}}{\sqrt{1-{\mathrm{\&rho;}}^2{\sin }^2}\mathrm{\&phi;}}+\cos \left(2\mathrm{\&alpha;}\right)\cos \mathrm{\&phi;}]$

The equal dimensionless of each variable in formula, the span of three parameters is defined as-0.9≤sin φ≤0.9,1,0≤ρ≤2.0 and-0.7≤sin α≤0.7;

To shot point x _swith imaging point (x ₀, T), make V _rmsfor the migration velocity at imaging point place, calculate

$\sin \mathrm{\&phi;}=\frac{x_0-x_s}{\sqrt{{(x_0-x_s)}^2+T^2{V}_{\mathrm{rms}}^2}}$

X in formula ₀and x _sunit be rice, when T is outward journey, unit second, V _rmsunit be meter per second.By rounding sin φ, imaging point place ρ and sin α to determining deviation of three the dimension directions shown, in table, obtain a according to the integer obtained ₁and a ₂; Make f _dbe the dominant frequency of seismic data, unit hertz, calculates

$c^{\&PlusMinus;}=-a_1\&PlusMinus;\sqrt{a_1^2-a_2/\left(2T{f}_d\right)}$

And then calculate two coordinate figures:

$x_r^{\&PlusMinus;}=x+c^{\&PlusMinus;}{\mathrm{TV}}_{\mathrm{rms}}\cos \mathrm{\&alpha;}-{\mathrm{TV}}_{\mathrm{rms}}\frac{\cos \mathrm{\&beta;}(\sin \mathrm{\&beta;}-p_3{c}^{\&PlusMinus;})(1-c^{\&PlusMinus;}\sin \mathrm{\&alpha;}/\mathrm{\&rho;})}{{\cos }^2\mathrm{\&beta;}+p_3\sin \mathrm{\&beta;}{c}^{\&PlusMinus;}-0.5p_3^2{\left(c^{\&PlusMinus;}\right)}^2}$

C in formula ^±dimensionless, unit be rice; Be calculated as follows shot point x again _s, imaging point (x ₀, T-1/ (2f _d)) and acceptance point corresponding incidence and reflection wave angle:

$\cos 2{\mathrm{\&phi;}}_{\mathrm{1,2}}=\frac{1}{\sqrt{1+p_s}\sqrt{1+p_r^{\&PlusMinus;}}}(1+0.5(p_s+p_r^{\&PlusMinus;}-\frac{{(x_s-x_r^{\&PlusMinus;})}^2}{{(T-1/\left(2f_d\right))}^2{V}_{\mathrm{rms}}^2}))$

In formula

$p_s=\frac{{(x_0-x_s)}^2}{V_{\mathrm{rms}}^2{T}^2},$ $p_r^{\&PlusMinus;}=\frac{{(x_0-x_r^{\&PlusMinus;})}^2}{V_{\mathrm{rms}}^2{T}^2}$

Just obtain imaging point (x ₀, T) and corresponding angle domain Fresnel zone [φ ₁, φ ₂], its unit degree of being.

Described does balancing energy to the two-dimensional migration road collection of each CDP point, first determine different time degree of depth upper angle territory Fresnel zone, then on angular deflection road collection, Fresnel zone part is excised, the remainder on Fresnel zone both sides is superposed respectively, obtain two imaging roads at each CDP point to be achieved in that and to do automatic gain to two-dimensional migration Dao Jizhu road, realize balancing energy; If the horizontal coordinate of CDP point is x ₀, the shot point coordinate of road centralized shotpoint dimension correspondence is x _s, time depth is T, α is the angle this CDP point different time degree of depth offseting lineups and surface level, take counter clockwise direction as positive-angle, V _rms(T) be the migration velocity in this CDP point different time degree of depth, Δ T is the sampling interval of time depth, calculates

$\sin \mathrm{\&phi;}=\frac{x_0-x_s}{\sqrt{{(x_0-x_s)}^2+T^2{V}_{\mathrm{rms}}^2\left(T\right)}},$ $\mathrm{\&rho;}=\frac{1}{\sqrt{\mathrm{\&Delta;T}}}\sqrt{T-(T-\mathrm{\&Delta;T})\frac{V_{\mathrm{rms}}^2(T-\mathrm{\&Delta;T})}{V_{\mathrm{rms}}^2\left(T\right)}}$

By setting up three-dimensional table spacing used in step 8, sin φ, ρ and sin α being rounded, in table, picking up two real numbers according to the integer obtained, and then calculate the angle domain Fresnel zone [φ in the different time degree of depth by the method for step 8 ₁, φ ₂], its unit degree of being;

In the different time degree of depth, the weight coefficient defined with angle change is

W(θ)＝1.0,θ≤φ ₁-5

$W\left(\mathrm{\&theta;}\right)=0.2+0.8\cos \left(\frac{\mathrm{\&pi;}(\mathrm{\&theta;}-{\mathrm{\&phi;}}_1)}{10}\right),$ φ ₁-5<θ<φ ₁

W(θ)=0.0,φ ₁≤θ≤φ ₂

$W\left(\mathrm{\&theta;}\right)=0.2+0.8\cos \left(\frac{\mathrm{\&pi;}(\mathrm{\&theta;}-{\mathrm{\&phi;}}_2)}{10}\right),$ φ ₂＜θ<φ ₂+5

W(θ)＝1.0,θ≥φ ₂+5

In above formula, θ is that skew road concentrates angle to tie up corresponding angle, unit degree.The angular-trace gather that each shot point is corresponding is concentrated to two-dimensional migration road, each time depth is multiplied by weight coefficient W (θ), just achieve excision Fresnel zone part.Angle corresponding to whole shot point is concentrated in two-dimensional migration road to be greater than and to be less than (φ ₁+ φ ₁the road of)/2 superposes respectively, just obtains two imaging roads at each CDP point.

Described on the basis at the lineups inclination angle that step 6 obtains, by the scan method that number percent increases and decreases, whether reflection line-ups is eliminated completely according on scattering wave migration stack section, determine to offset accurately lineups inclination angle and be achieved in that the angle of lineups and surface level on the migration stack section that obtain step 6 carries out number percent increase and decrease, to the universe angle value applying step 9 that each number percent increases and decreases, and by two imaging trace-stackings that each CDP point place obtains, form stacked section; The migration stack Profile Correlation this section and step 5 obtained, if disappear at a certain imaging point reflection line-ups, angular values is now exactly these lineups inclination angles accurately.

The described accurate skew lineups inclination angle utilizing step 10 to obtain, repeat step 9, obtain two migration stack sections be made up of imaging road on the right side of imaging road on the left of Fresnel zone and Fresnel zone respectively, this two section applied energies equilibrium is subtracted each other, obtain only being achieved in that the scattering wave migration stack section of the independent imaging of little turn-off breakpoint and utilize reflecting interface and horizontal plane angle accurately, namely lineups inclination angle is offset, repeat step 9, the imaging road on the left of each CDP point Fresnel zone and right side is combined respectively, two the migration imaging sections obtained, remember that its amplitude is A (x respectively, and B (x T), T), wherein x is the horizontal coordinate of CDP point, T is time depth, right | A (x, T) | with | B (x, T) | carry out large scale smoothly, obtain its low frequency component, be designated as E (x, T) and F (x, T), to the scattering wave migration stack section of the independent imaging of little turn-off breakpoint be then only

$I(x,T)=(\frac{A(x,T)}{E(x,T)+\mathrm{\&epsiv;}}-\frac{B(x,T)}{F(x,T)+\mathrm{\&epsiv;}})\sqrt{E(x,T)F(x,T)}$

In above formula, ε is a stable factor a small amount of.

The scattering wave Prestack Imaging method of the little turn-off breakpoint of identification of the present invention, can obtain the distributed image of the little turn-off breakpoint of subsurface geological structure.

The scattering wave Prestack Imaging method of the little turn-off breakpoint of identification of the present invention, can be less than the independent imaging of breakpoint of 1/4 wavelength to the impalpable turn-off of typical reflection ripple imaging, can significantly improve the recognition capability of method of seismic prospecting to crack, underground and tomography.

The scattering wave Prestack Imaging method of the little turn-off breakpoint of identification of the present invention, can suppress transmitted wave energy in migration process, outstanding breakpoint scattering wave.

The scattering wave Prestack Imaging method of the little turn-off breakpoint of identification of the present invention, independently can determine the migration velocity that scattering wave skew uses and interface dip.

Specific implementation principle of the present invention is as follows:

The present invention adopts and concentrates depressor reflex ripple to realize scattering wave imaging in migration imaging road, and utilizes the reversal of poles of little breakpoint scattering wave, strengthens little breakpoint imaging; Its core offsets the different manifestations on road collection according to reflection line-ups and scattering wave lineups at single big gun, carrys out depressor reflex ripple, retains scattering wave.Mainly contain following 3 points, one is improve prestack time migration method, obtains the angle domain skew road collection of big gun record; Two is the Fresnel zones accurately determining reflection wave on angle domain skew road collection, takes this depressor reflex ripple; Three is breakpoint scattering wave image enhancements.

Acquisition and processing carries out in two dimensional cross-section, therefore will carry out discussion for two-dimensional problems.Specific implementation principle is as follows:

1. angle domain skew road collection

Prestack time can regard that pre-stack depth migration is in velocity equivalent parameter as, being similar to namely under root-mean-square velocity.Based on phase-shift method and the steady phase point principle of pre-stack depth migration, conventional pre-stack time migration formula can be derived, and based on this framework, can when not needing known interface normal, obtain the angle of incident wave and interface normal, the angle domain taking this to obtain incident wave corresponding to single shot record and interface normal angle offsets road collection.

In wave number-frequency field, each seismic trace of big gun record can be expressed as F _m(ω) exp (-jk _xx _g), wherein ω is angular frequency, k _xhorizontal wave number, x _gbe the horizontal coordinate of this seismic trace acceptance point, j is unit imaginary number, F _m(ω) be the Fourier transform of this seismic trace time-domain signal.

If single seismic trace to be regarded as a big gun record, further nonhomogeneous media is approximately layered medium, by the degree of depth-wave number-frequency field phase-shift method can the degree of depth continuation result of this seismic trace be

$\tilde{P}(k_x,\mathrm{\&omega;},T=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\mathrm{\&Delta;}{T}_i)=F_m\left(\mathrm{\&omega;}\right)\exp (-j{k}_x{x}_g)\exp \left(\mathrm{j\&omega;}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\mathrm{\&Delta;}{T}_i\sqrt{1-\frac{v_i^2}{{\mathrm{\&omega;}}^2}{k}_x^2}\right)---\left(1\right)$

The degree of depth is expressed with the straight whilst on tour T of one way that hangs down in formula, Δ T _ithe thickness of expressing when being each layer medium outward journey, has Δ z _i=v _iΔ T _i, v _ibe the speed of each layer medium, n is the medium number of plies of more than zone of interest, it is the wave number-frequency field continuation result at time depth T place.

Increment part in the exponential term of formula (1) Exponential function can approximate expression be:

$\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left(\mathrm{\&Delta;}{T}_i\sqrt{1-k_x^2{v}_i^2/{\mathrm{\&omega;}}^2}\right)\&ap;\sqrt{1-k_x^2{V}_{\mathrm{rms}}^2/{\mathrm{\&omega;}}^2}\&CenterDot;\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\mathrm{\&Delta;}{T}_i\right)---\left(2\right)$

Be Taylor to (2) formula both sides to launch, and be truncated to Section 2, can obtain:

$V_{\mathrm{rms}}^2=\frac{1}{T}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{v}_i^2\mathrm{\&Delta;}{T}_i---\left(3\right)$

V _rmsnamely be the migration velocity required for pre-stack time migration, namely conventional root-mean-square velocity.

Formula (2) is substituted into (1) formula and does space inverse fourier transform, the space-frequency field continuation result obtaining single seismic trace is:

$P(x,\mathrm{\&omega;},T)=\frac{{\mathrm{\&omega;}}^2}{4{\mathrm{\&pi;}}^2}\&Integral;F_m\left(\mathrm{\&omega;}\right)\exp \left\{\mathrm{j\&omega;}(T\sqrt{1-V_{\mathrm{rms}}^2{p}_x^2}+p_x(x-x_g))\right\}{\mathrm{dp}}_x---\left(4\right)$

Ray parameter p is introduced in formula _x=k _x/ ω.The available steady phase point principle of formula (4) tries to achieve progressive solution, and the ray parameter that surely phase point is corresponding then represents the wavefront surface normal direction of anti-pass wave field corresponding to continuation result.The steady phase point of formula (4) can be solved by following formula:

$\frac{\mathrm{d\&phi;}\left(p_x\right)}{{\mathrm{dp}}_x}=0---\left(5\right)$

Formula has in (5)

$\mathrm{\&phi;}\left(p_x\right)=T\sqrt{1-V_{\mathrm{rms}}^2{p}_x^2}+p_x(x-x_g)$

Can be solved by formula (5):

$p_x=\frac{x-x_g}{V_{\mathrm{rms}}\sqrt{{\left({x-x}_g\right)}^2+V_{\mathrm{rms}}^2{T}^2}}---\left(6\right)$

Formula (6) is substituted into steady phase point principle formula, and the integration that can obtain formula (4) is

$P(x,\mathrm{\&omega;},T)=F_m\left(\mathrm{\&omega;}\right)\sqrt{\frac{\mathrm{\&omega;}}{2\mathrm{\&pi;}}}\exp (-j\frac{\mathrm{\&pi;}}{4})\frac{T}{{\mathrm{\&tau;}}_g^{3/2}{V}_{\mathrm{rms}}}\exp \left(\mathrm{j\&omega;}{\mathrm{\&tau;}}_g\right)---\left(7\right)$

τ in formula _gfor seimic travel time, have

${\mathrm{\&tau;}}_g=T\sqrt{1+\frac{{(x-x_g)}^2}{V_{\mathrm{rms}}^2{T}^2}}---\left(8\right)$

Identical with conventional pre-stack time migration during the walking of formula (8).Because actual focus is point source, the geometrical attenuation item of seismic event need consider seismic event propagation in three dimensions, and the amplitude item therefore in formula (7) should utilize three-dimensional amplitude formula, and this pattern (7) can be rewritten as

$P(x,\mathrm{\&omega;},T)=F_m\left(\mathrm{\&omega;}\right)\sqrt{\frac{\mathrm{\&omega;}}{2\mathrm{\&pi;}}}\exp (-j\frac{\mathrm{\&pi;}}{4})\frac{T}{{\mathrm{\&tau;}}_g^2{V}_{\mathrm{rms}}^2}\exp \left(j{\mathrm{\&omega;\&tau;}}_g\right)---\left(9\right)$

Adopt same procedure, the ray parameter that can obtain the shot point incident wave of this single shot record is

$p_x=\frac{x-x_s}{V_{\mathrm{rms}}\sqrt{{(x-x_s)}^2+V_{\mathrm{rms}}^2{T}^2}}---\left(10\right)$

And the main story wave field of shot point incident wave is in the different time degree of depth

$P_S(x,y,\mathrm{\&omega;},T)=S\left(\mathrm{\&omega;}\right)\sqrt{\frac{\mathrm{\&omega;}}{2\mathrm{\&pi;}}}\exp \left(j\frac{\mathrm{\&pi;}}{4}\right)\frac{T}{{\mathrm{\&tau;}}_s^2{V}_{\mathrm{rms}}^2}\exp (-j{\mathrm{\&omega;\&tau;}}_s)---\left(11\right)$

In formula, S (ω) is the Fourier transform of source wavelet, τ _sfor focus is to the seimic travel time of imaging point, have

${\mathrm{\&tau;}}_s=T\sqrt{1+\frac{{(x-x_s)}^2}{V_{\mathrm{rms}}^2{T}^2}}---\left(12\right)$

Generally source wavelet is unknown.Since deconvolution can reject the impact of wavelet, we can ignore the impact of source wavelet continuous item in imaging, namely if establish big gun to record n seismic trace, application deconvolution image-forming condition, the imaging results that can obtain single shot record is

$I(x,T)=\underset{m=1}{\overset{n}{\mathrm{\&Sigma;}}}{({\mathrm{\&tau;}}_s/{\mathrm{\&tau;}}_g^m)}^2\&Integral;F_m\left(\mathrm{\&omega;}\right)\sqrt{\frac{\mathrm{\&omega;}}{2\mathrm{\&pi;}}}\exp (-j\frac{\mathrm{\&pi;}}{4})\exp \left(\mathrm{j\&omega;}({\mathrm{\&tau;}}_s+{\mathrm{\&tau;}}_g^m)\right)\mathrm{d\&omega;}---\left(13\right)$

$=\underset{m=1}{\overset{n}{\mathrm{\&Sigma;}}}{({\mathrm{\&tau;}}_s/{\mathrm{\&tau;}}_g^m)}^2{g}_m({\mathrm{\&tau;}}_s+{\mathrm{\&tau;}}_g^m)$

In formula it is half derivative of seismic trace time-domain signal; be the acceptance point of m seismic trace to the walking of imaging point time, obtain by the horizontal coordinate of this acceptance point is substituted into formula (8).

If the interval velocity assuming picture point (x, T) place medium is v, according to the ray parameter of formula (6) and (10), more than the wavefront surface direction that can obtain shot point incident field and acceptance point anti-pass wave field, device for carrying a tripot is vp _xwith and then shot point incident wave corresponding to this seismic trace can be obtained and acceptance point anti-pass ripple at the angle of imaging point is

$\cos 2\mathrm{\&theta;}=\frac{1}{\sqrt{1+p_s}\sqrt{1+p_r}}(\sqrt{1+(1-{\mathrm{\&rho;}}^2)p_s}\sqrt{1+(1-{\mathrm{\&rho;}}^2)p_r}+0.5(p_s+p_r-\frac{H^2}{T^2{V}_{\mathrm{rms}}^2}))---\left(14\right)$

In formula, 2 θ are angles of incident wave and reflection wave, ρ=v/V _rms, H is the offset distance of seismic trace, parameter p _sand p _rfor

$p_s=\frac{{(x-x_s)}^2}{V_{\mathrm{rms}}^2{T}^2},$ $p_r=\frac{{(x-x_g)}^2}{V_{\mathrm{rms}}^2{T}^2}$

For distinguishing reflection and scattering wave better, the present invention introduces positive and negative to the angle that formula (14) is tried to achieve, definition x _g>=x _stime angle be just.

Practical application middle level speed v is unknown.Because angle gathers of the present invention is only for depressor reflex ripple, instead of utilizing angle-magnitude relation to carry out inverting, therefore can suppose ρ=1 when determining angle θ, this pattern (14) can approximate expression be

$\cos 2\mathrm{\&theta;}=\frac{1}{\sqrt{1+p_s}\sqrt{1+p_r}}(1+0.5(p_s+p_r-\frac{H^2}{T^2{V}_{\mathrm{rms}}^2}))---\left(15\right)$

If definition Δ θ is angle varied pitch, formula (13) cumulative in each imaging point (x, T) angle that each seismic trace is corresponding is considered, by k Δ θ-Δ θ/2≤θ <k Δ θ+Δ θ/2, each angular range is added up respectively, just can obtain the angle domain skew road collection I (x of single shot record, T, k Δ θ).For reducing the calculated amount directly being calculated angle by formula (15), cos can be obtained by formula (15) ²θ, by itself and cos ²(k Δ θ ± Δ θ/2) compare, and determine the cumulative position of this skew amplitude.Work as x _r>=x _s, skew amplitude is added in positive-angle, otherwise, be then added to negative angle.

2. angle domain Fresnel zone

Adding up along angle to angle domain skew road collection I (x, T, θ), obtaining in the process of the migration result I (x, T) of formula (13), reflecting the contribution of wave-wave relative to imaging point mainly from the Fresnel zone at steady phase point place; And the angle that surely phase point is corresponding is exactly the angle of shot point incident wave and reflecting interface normal.First excise Fresnel zone part if can offset to concentrate in angle domain and then superpose, the imaging of depressor reflex ripple can be realized.Therefore obtaining Fresnel zone is accurately a key of the present invention.

So-called Fresnel zone, be namely in this region, the time depth difference that the skew lineups of reflection wave and imaging point place offset lineups is less than specific spacing, such as 1/4 cycle; Because seismic signal is in a frequency band, this time interval is defined as 1/2 cycle corresponding to seismic event dominant frequency by the present invention.Based on the ray parameter expression formula of formula (6) and (10), the angle that the Fresnel zone two ends at inclined reflection interface are corresponding on angle gathers will be tried to achieve, i.e. angle domain Fresnel zone below.

As shown in Figure 1, if the inclination angle of reflecting interface is α, the coordinate of shot point S is x _s, the coordinate of imaging point I is (x, T), and it is also the incidence point of incident wave on reflecting interface; If on interface, the point being d apart from imaging point along interface is the border of Fresnel zone, i.e. I ' point, then the picture of this point reflection ripple in imaging point angle gathers will at (x, T-δ) place, i.e. P point, wherein δ=1/ (2f _d), f _dit is the dominant frequency of seismic data; If I ' point reflection ripple shines earth's surface R point, its coordinate is x _rpoint, by incident wave identical with the angle of interface normal with reflection wave (intersecting at I ' point), and shot point arrives x again to Fresnel zone frontier point (I ' point) _rthese two conditions identical when arriving walking of R point with shot point again to P point during the walking of point, can obtain following two equations:

$\frac{1}{\sqrt{{(x+d\cos \mathrm{\&alpha;}-x_s)}^2+{(V_{\mathrm{rms}}T-d\sin \mathrm{\&alpha;}/\mathrm{\&rho;})}^2}}\left(\mathrm{\&rho;}\sin \mathrm{\&alpha;}(x+d\cos \mathrm{\&alpha;}-x_s\right)+$

$\cos \mathrm{\&alpha;}\sqrt{{(V_{\mathrm{rms}}T-d\sin \mathrm{\&alpha;}/\mathrm{\&rho;})}^2+(1-{\mathrm{\&rho;}}^2){(x+d\cos \mathrm{\&alpha;}-x_s)}^2})=\frac{1}{\sqrt{{(x-d\cos \mathrm{\&alpha;}-x_r)}^2+{(V_{\mathrm{rms}}T-d\sin \mathrm{\&alpha;}/\mathrm{\&rho;})}^2}}$

$(\mathrm{\&rho;}\sin \mathrm{\&alpha;}(x+d\cos \mathrm{\&alpha;}-x_r)+\cos \mathrm{\&alpha;}\sqrt{{(V_{\mathrm{rms}}T-d\sin \mathrm{\&alpha;}/\mathrm{\&rho;})}^2+(1-{\mathrm{\&rho;}}^2){(x-d\cos \mathrm{\&alpha;}-x_r)}^2})---\left(16\right)$

$(T+d\sin \mathrm{\&alpha;}/\left(\mathrm{\&rho;}{V}_{\mathrm{rms}}\right))(\sqrt{1+\frac{{(x+d\cos \mathrm{\&alpha;}-x_s)}^2}{{(V_{\mathrm{rms}}T-d\sin \mathrm{\&alpha;}/\mathrm{\&rho;})}^2}}+\sqrt{1+\frac{{(x+d\cos \mathrm{\&alpha;}-x_r)}^2}{{(V_{\mathrm{rms}}T-d\sin \mathrm{\&alpha;}/\mathrm{\&rho;})}^2}})---\left(17\right)$

$=(T-1/\left(2f_d\right))(\sqrt{1+\frac{{(x-x_s)}^2}{V_{\mathrm{rms}}^2{(T-1/\left(2f_d\right))}^2}}+\sqrt{1+\frac{{(x-x_r)}^2}{V_{\mathrm{rms}}^2{(T-1/\left(2f_d\right))}^2}})$

ρ=v/V in formula _rms, d is the distance along imaging point Fresnel zone border, direction, interface, x _rit is the acceptance point coordinate intending solving.Two groups of (x are solved by formula (16) and (17) _r, d), by x _s, x _r(x, T-1/ (2f _d)) substitute into the angle that formula (15) can try to achieve Fresnel zone two ends correspondence on angle gathers.

In formula (16) and (17), ρ=1 can not be supposed.And v can be solved by formula (3) and be

$v^2=\frac{1}{\mathrm{\&Delta;T}}[{\mathrm{TV}}_{\mathrm{rms}}^2(x,T)-(T-\mathrm{\&Delta;T})V_{\mathrm{rms}}^2(x,T-\mathrm{\&Delta;T})]$

Thus

${\mathrm{\&rho;}}^2=\frac{1}{\mathrm{\&Delta;T}}[T-(T-\mathrm{\&Delta;T})\frac{V_{\mathrm{rms}}^2(x,T-\mathrm{\&Delta;T})}{V_{\mathrm{rms}}^2(x,T)}]---\left(18\right)$

Sampling interval Δ T due to time depth in practical application is less, and for avoiding error, the Δ T in formula (18) can substitute with 2 Δ T or 3 Δ T.When calculating formula (18), should first to V _rmsfully level and smooth.

Direct solution formula (16) and (17) this group nonlinear equation are difficult, utilize d and TV in formula (16) _rmscomparing is in a small amount, and 1/ (2f in formula (17) _d) be in a small amount with T-phase ratio, by ignoring d/ (TV _rms) and 1/ (2Tf _d) higher order term, simplified style (16) and (17), obtain the approximate solution of degree of precision.

Actual computation can adopt look-up table to try to achieve x _rtwo solutions.Set up three-dimensional table with sin φ, ρ and sin α for dimension, deposit by two real numbers of trying to achieve as shown in the formula (19) and (20) in table:

$a_1=\frac{\mathrm{\&rho;}(p_1\sin \mathrm{\&phi;}\cos \mathrm{\&beta;}+p_2\sin \mathrm{\&beta;}\cos \mathrm{\&phi;})-\sin \mathrm{\&alpha;}(\cos \mathrm{\&beta;}+\cos \mathrm{\&phi;})}{\mathrm{\&rho;}\cos \mathrm{\&phi;}\cos \mathrm{\&beta;}(p_1^2\cos \mathrm{\&phi;}-p_2^2\cos \mathrm{\&beta;})-2p_1\sin \mathrm{\&alpha;}\sin \mathrm{\&phi;}\cos \mathrm{\&beta;}+\left[2p_3\cos \mathrm{\&phi;}(\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}-\mathrm{\&rho;}{p}_2)\right]/{\cos }^2\mathrm{\&beta;}}---\left(19\right)$

$a_2=\frac{\mathrm{\&rho;}\cos \mathrm{\&phi;}\cos \mathrm{\&beta;}(\cos \mathrm{\&phi;}+\cos \mathrm{\&beta;})}{\mathrm{\&rho;}\cos \mathrm{\&phi;}\cos \mathrm{\&beta;}(p_1^2\cos \mathrm{\&phi;}-p_2^2\cos \mathrm{\&beta;})-2p_1\sin \mathrm{\&alpha;}\sin \mathrm{\&phi;}\cos \mathrm{\&beta;}+\left[2p_3\cos \mathrm{\&phi;}(\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}-\mathrm{\&rho;}{p}_2)\right]/{\cos }^2\mathrm{\&beta;}}---\left(20\right)$

In upper two formulas

$\sin \mathrm{\&beta;}=\sin \left(2\mathrm{\&alpha;}\right)\sqrt{1/p^2-{\sin }^2\mathrm{\&phi;}}-\cos \left(2\mathrm{\&alpha;}\right)\sin \mathrm{\&phi;},$ $\cos \mathrm{\&beta;}=\sqrt{1-{\sin }^2\mathrm{\&beta;}},$ $\cos \mathrm{\&phi;}=\sqrt{1-{\sin }^2\mathrm{\&phi;}}$

p ₁＝cosαcosφ+sinαsinφ/ρ，p ₂＝cosαcosβ+sinαsinβ/ρ

$p_3=p_1\cos \mathrm{\&phi;}[\frac{\mathrm{\&rho;}\sin \left(2\mathrm{\&alpha;}\right)\sin \mathrm{\&phi;}\cos \mathrm{\&phi;}}{\sqrt{1-{\mathrm{\&rho;}}^2{\sin }^2}\mathrm{\&phi;}}+\cos \left(2\mathrm{\&alpha;}\right)\cos \mathrm{\&phi;}]$

To arbitrary shot point and imaging point, calculate,

$\sin \mathrm{\&phi;}=\frac{x-x_s}{\sqrt{{(x-x_s)}^2+T^2{V}_{\mathrm{rms}}^2}}$

According to the ρ at sin φ, imaging point place and sin α corresponding to inclination alpha, in table, pick up a ₁and a ₂; Calculate

$c^{\&PlusMinus;}=-a_1\&PlusMinus;\sqrt{a_1^2-a_2/\left(2T{f}_d\right)}$

X can be obtained _rtwo solutions

$x_r^{\&PlusMinus;}=x+c^{\&PlusMinus;}{\mathrm{TV}}_{\mathrm{rms}}\cos \mathrm{\&alpha;}-{\mathrm{TV}}_{\mathrm{rms}}\frac{\cos \mathrm{\&beta;}(\sin \mathrm{\&beta;}-p_3{c}^{\&PlusMinus;})(1-c^{\&PlusMinus;}\sin \mathrm{\&alpha;}/\mathrm{\&rho;})}{{\cos }^2\mathrm{\&beta;}+p_3\sin \mathrm{\&beta;}{c}^{\&PlusMinus;}-0.5p_3^2{\left(c^{\&PlusMinus;}\right)}^2}---\left(21\right)$

By formula (21) and imaging point (x, T-1/ (2f _d)) substitute into formula (15), two angle φ that Fresnel zone is corresponding can be solved ₁and φ ₂.

If excision angular range [φ in angle domain skew road collection I (x, T, θ) ₁, φ ₂] in value, and then along angle add up, then will not there is reflection line-ups in migration result I (x, T).This is because the value in the Fresnel zone contributed to some extent reflection wave is cut.Excise and realize by introducing weight coefficient in cumulative process, for avoiding edge effect, weight coefficient should be transitioned into 0 by 1 glossily.

When determining Fresnel zone, the inclination angle of reflecting interface need be utilized, i.e. the angle of reflection line-ups and surface level; What obtain due to prestack time migration method is time domain migrated section, generally only can obtain approximate reflector dip based on this section.Fresnel zone scope is vital to excision reflection wave accurately, and from formula (19)-(21), interface dip has material impact to Fresnel zone scope.For addressing this problem, the scan method of the present invention's application number percent increase and decrease obtains accurate reflector dip.The foundation of scanning is exactly that can scattering wave imaging section reject reflection line-ups completely.Whether this, by contrast migration stack section, observes continuous print reflection line-ups and disappears to realize.

3. breakpoint scattering wave image enhancement

Added up along angle by the weight coefficient introduced for Fresnel zone, in theory, can obtain not containing single big gun migration result of reflection wave; All single big gun migration result add up, and can obtain not containing the scattering wave migration stack section of reflection wave.But when scattering wave is more weak or the reception number of channels of single shot record is less, reflection wave remaining after adding up along angle and noise be magnitude same with scattering wave amplitude likely, be so just difficult to obtain scattering wave imaging results clearly.For this reason, need strengthen scattering wave.By doing automatic gain to the skew Dao Jizhu road before adding up along angle, realizing balancing energy, can strengthen scattering wave on each seismic trace, because scattering wave has good same phasic property, superposition road concentrates scattering wave to be enhanced, and noise can't increase.Such process can strengthen the energy of scattering wave in migration stack section.The angle of the present invention to incident wave and reflection wave is introduced positive and negative, by extremely efficient angle, improves the actual effect of this enhancing process.

The present invention carries out imaging by utilizing scattering wave to little breakpoint, thus the minor fault of better understanding underground structure and breakpoint.Be different from medium small scale and change the scattering wave caused, a key character of breakpoint scattering wave is scattering wave lineups by there is reversal of poles during steady phase point.Utilizing this feature, by subtracting each other the imaging results of steady phase point both sides, giving prominence to reflection wave further; This subtraction can suppress the scattering wave imaging of other type and remaining reflection line-ups, plays an important role to the little breakpoint imaging of lifting.Its specific implementation process is as follows:

Concentrate the road of Fresnel zone both sides to superpose respectively in two-dimensional migration road, obtain two imaging roads at each CDP point.And then two migration stack sections be made up of imaging road on the right side of imaging road on the left of Fresnel zone and Fresnel zone respectively can be obtained.The amplitude of note two migration stack sections is A (x, T) and B (x, T) respectively, and wherein x is the horizontal coordinate of CDP point, and T is time depth.Right | A (x, T) | with | B (x, T) | carry out large scale smoothly, obtain its low frequency component, be designated as E (x, T) and F (x, T).To the scattering wave migration stack section of the independent imaging of little turn-off breakpoint be then only

$I(x,T)=(\frac{A(x,T)}{E(x,T)+\mathrm{\&epsiv;}}-\frac{B(x,T)}{F(x,T)+\mathrm{\&epsiv;}})\sqrt{E(x,T)F(x,T)}---\left(22\right)$

In above formula, ε is a stable factor a small amount of.

Beneficial effect of the present invention: the method uses conventional reflected seismic information, can to the independent imaging of the impalpable little turn-off breakpoint of typical reflection ripple offset method, can significantly improve the recognition capability of method of seismic prospecting to crack, underground and tomography.The method can indicate minor fault, the breakpoint of underground structure subtly, thus can be familiar with migration and the shutoff of underground oil and gas better, has significant application value to the exploration and development of hydrocarbon resources.

Accompanying drawing explanation

Fig. 1 is the schematic diagram determining Fresnel zone.Figure middle conductor AB represents the reflecting interface of inclination, and α is the level inclination at interface, and S is shot point, and imaging point is I (x, T), is incident to the reflection wave outgoing of I ' to surperficial R point, earthquake wave trajectory in represented by dotted arrows nonhomogeneous media; If the seismic event that R point receives can in P (x, T-δ) imaging by principle when waiting, and δ was 1/2 cycle, then namely I ' be the Fresnel zone frontier point of imaging point I.

Fig. 2 is the isogram of the South Sea block migration velocity field obtained by actual seismic data.In figure, numeral is migration velocity value.In figure, horizontal coordinate presses the overall coordinate definition of block, therefore does not start from scratch.

Fig. 3 is the migration stack section that block conventional pre-stack time migration in the South Sea obtains.

Fig. 4 is the angle dimension tangent plane of the two-dimensional migration road collection of the corresponding shot point coordinate 16.45km of CDP point 700.Visible reflectance ripple lineups show bending shape.Middle blank causes owing to lacking near migration range seismologic record.

Fig. 5 is the result of Fig. 4 angular deflection road collection excision Fresnel zone.

Fig. 6 a and Fig. 6 b is that before and after Fig. 4 angular deflection road collection excision Fresnel zone, stack result compares.

Fig. 6 a is the stack result before excision Fresnel zone, and Fig. 6 b is the stack result after excision Fresnel zone, and visible original stronger reflex amplitude is pressed.

Fig. 7 is the migration stack section that on the left of Fresnel zone, imaging road is formed.The conventional migration technique stacked section of comparison diagram 3, visible reflectance ripple lineups are substantially disallowable.

Fig. 8 is the migration stack section that on the right side of Fresnel zone, imaging road is formed.The conventional migration technique stacked section of comparison diagram 3, visible reflectance ripple lineups are substantially disallowable.

Fig. 9 is only to the scattering wave migration stack section of the independent imaging of little turn-off breakpoint.Comparison diagram 7 and Fig. 8 visible, remaining reflection line-ups and other interference are suppressed completely, and the breakpoint of tomography is strengthened.

Figure 10 a and Figure 10 b gives the little breakpoint distributed image of Fig. 9 and Fig. 3 conventional migration technique stacked section at horizontal coordinate 13.2km, the local contrast of time depth 1.0s neighborhood.Figure 10 a is the partial enlargement of Fig. 3 conventional migration technique result, and Figure 10 b is the partial enlargement of the little breakpoint distribution plan of Fig. 9, contrasts, occur the position of little change at reflection line-ups from figure, and the imaging of breakpoint scattering wave just indicates breakpoint to be existed.

Figure 11 a and Figure 11 b gives the little breakpoint distributed image of Fig. 9 and Fig. 3 conventional migration technique stacked section at horizontal coordinate 13.2km, the local contrast of time depth 1.5s neighborhood.Figure 11 a is the partial enlargement of Fig. 3 conventional migration technique result, and Figure 11 b is the partial enlargement of the little breakpoint distribution plan of Fig. 9; The same with Figure 10, occur the position of little change at reflection line-ups, the imaging of breakpoint scattering wave just indicates the existence of breakpoint.

Figure 12 gives the little breakpoint distributed image of Fig. 9 and Fig. 3 conventional migration technique stacked section at horizontal coordinate 14.7km, the local contrast of time depth 2.4s neighborhood.Figure 12 a is the partial enlargement of Fig. 3 conventional migration technique result, and Figure 12 b is the partial enlargement of the little breakpoint distribution plan of Fig. 9.In Figure 12 a, the change of reflection line-ups is not obvious, but the imaging of Figure 12 b breakpoint scattering wave explicitly points out the existence of breakpoint.This shows the impalpable little turn-off breakpoint of this anti-bright identifiable design typical reflection ripple skew.

Embodiment

Embodiment 1: the scattering wave Prestack Imaging method identifying little turn-off breakpoint, is example for South Sea block, is specially following steps:

(1) with the seismic reflection signals of the artificial epicenter excitation of acquisition mode record of single towing cable, be recorded on tape.Concrete seismic reflection signals is, adopts wall scroll towing cable collection, every bar towing cable 944 road, track pitch 3.125m; In the front end of towing cable with air gun battle array epicenter excitation, the distance of air gun battle array and first reception channel is 125m, record length 4s, time-sampling 1ms; Excite altogether and record 680 big guns, shot point spacing 12.5m.

(2) read seismic signal from tape, form Prestack seismic data, conventional sound attenuation process is done to Prestack seismic data; Build a two dimensional cross-section by along line direction and depth direction, two dimensional processing methods is used to seismic data; For part common midpoint, extract common midpoint gather, conventional NMO velocity pickup is done to the common midpoint gather extracted, the reciprocal value of obtained speed is done laterally average, using the laterally uniform velocity field that obtains as initial offset speed.

(3) according to initial offset speed, carry out calculations of offset by conventional prestack time method, form CRP gather.

(4) to CRP gather, utilize initial offset speed to do inverse dynamic correction, then do normal moveout correction and obtain the speed after upgrading, do space smoothing to the reciprocal value of speed after the renewal of each CDP point, namely the level and smooth velocity field obtained is migration velocity field.Fig. 2 is the isogram of migration velocity field.In figure, numeral is migration velocity value.

(5) according to migration velocity field, conventional pre-stack time migration is done to Prestack seismic data and calculates, the dynamic school of high-order residue and the excision that stretches are done to the CRP gather obtained, the start offset distance of record residual normal moveout and the excision that stretches; Along offset distance superposition CRP gather, obtain migration stack section.Fig. 3 is the migration stack section that conventional pre-stack time migration obtains.

(6) according to migration stack section and migration velocity field, determining the angle of lineups and surface level on migrated section, take counter clockwise direction as positive-angle.

(7) according to migration velocity field, big gun territory pre-stack time migration is done to Prestack seismic data, form the two-dimensional migration road collection with shot point and angle change at each CDP point.Fig. 4 is the skew road collection of the big gun record of CDP point 700 place shot point coordinate 16.45km, i.e. an angle dimension tangent plane of two-dimensional migration road collection.This big gun record is by the single shot record expanded as bilateral record.

To each CDP point, define a three-dimensional array and deposit migration result, the corresponding time depth of one dimension, another dimension is shot point numbering, and the third dimension is angle, and angle is from bearing 60 degree to positive 60 degree of changes; Prestack seismic data is pressed common-shot-gather sequence; To the common-shot-gather of monolateral record, utilize big gun, the cautious principle such as mutual of exchanging, the common-shot-gather of the bilateral record of constructing virtual; By shot point coordinate, common-shot-gather is divided into groups, difference group common-shot-gather is placed on the different computing nodes of cluster computer; At each computing node, whole seismic trace is circulated; To each seismic trace, to the whole CDP dot cycles in migration aperture; To each CDP point, determine initial imaging point by migration aperture, the imaging point below this imaging point is calculated successively; First by shot point, geophone station and imaging point coordinate and migration velocity field, when determining to walk and imaging weight coefficient, and then determined to offset amplitude A by the seismologic record amplitude of this seismic trace; Then according to the offset distance of seismic trace, the horizontal coordinate of imaging point and time depth, the residual normal moveout that pickup step 5 records and the start offset distance stretching excision; The time depth of foundation offset distance and imaging point, judge whether skew amplitude A should be cut, if do not excise, then the angle calculating shot point incident wave and reflecting interface normal is

${\cos }^2\mathrm{\&theta;}=\frac{0.5}{\sqrt{1+p_s}\sqrt{1+p_r}}(1+0.5(p_s+p_r-\frac{H^2}{T^2{V}_{\mathrm{rms}}^2}))+0.5$

In above formula, θ is the angle of shot point incident wave and reflecting interface normal, unit degree, and T is the time depth of imaging point, namely during outward journey, and unit second, V _rmsfor imaging point (x ₀, T) and the migration velocity at place, unit meter per second, H is the offset distance of seismic trace, unit rice, p _sand p _rfor

$p_s=\frac{{(x_0-x_s)}^2}{V_{\mathrm{rms}}^2{T}^2},$ $p_r=\frac{{(x_0-x_g)}^2}{V_{\mathrm{rms}}^2{T}^2}$

Wherein x ₀, x _gand x _sthe horizontal coordinate of imaging point, geophone station and shot point respectively, unit rice;

If Δ θ is the angular interval of predefined, unit degree, according to cos ²(k Δ θ-Δ θ 2)≤cos ²θ <cos ²(k Δ θ-Δ θ/2) determine positive integer k; If the big gun belonging to this seismic trace is numbered n, if the residual normal moveout of pickup is τ, its unit is second, works as x _g>=x _s, then skew amplitude A is added on (T+ τ, n, the k Δ θ) element of the three-dimensional array depositing migration result, namely in positive-angle; Work as x _g<x _s, then skew amplitude A is added on (T+ τ, the n ,-k Δ θ) element of the three-dimensional array depositing migration result; Complete to whole seismic trace circulation, just obtain the two-dimensional migration road collection with shot point and angle change at each CDP point, this angle is the angle of shot point incident wave and reflecting interface normal.

(8) according to the angle of migration velocity field, shot point coordinate, imaging point coordinate, imaging point place skew lineups and surface level, set up the three-dimensional table determining Fresnel zone, three-dimensional table determines the angle domain Fresnel zone of each imaging point accordingly.

Make ρ be the ratio of imaging point place actual layer speed and migration velocity, α is the inclination angle of imaging point place reflecting interface, and namely offseting the angle of lineups and surface level, take counter clockwise direction as positive-angle; Be that three dimensions set up three-dimensional table with sin φ, ρ and sin α, deposit two dimensionless factors in table, its value is calculated by following two formulas:

$a_1=\frac{\mathrm{\&rho;}(p_1\sin \mathrm{\&phi;}\cos \mathrm{\&beta;}+p_2\sin \mathrm{\&beta;}\cos \mathrm{\&phi;})-\sin \mathrm{\&alpha;}(\cos \mathrm{\&beta;}+\cos \mathrm{\&phi;})}{\mathrm{\&rho;}\cos \mathrm{\&phi;}\cos \mathrm{\&beta;}(p_1^2\cos \mathrm{\&phi;}-p_2^2\cos \mathrm{\&beta;})-2p_1\sin \mathrm{\&alpha;}\sin \mathrm{\&phi;}\cos \mathrm{\&beta;}+\left[2p_3\cos \mathrm{\&phi;}(\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}-\mathrm{\&rho;}{p}_2)\right]/{\cos }^2\mathrm{\&beta;}}$

$a_2=\frac{\mathrm{\&rho;}\cos \mathrm{\&phi;}\cos \mathrm{\&beta;}(\cos \mathrm{\&phi;}+\cos \mathrm{\&beta;})}{\mathrm{\&rho;}\cos \mathrm{\&phi;}\cos \mathrm{\&beta;}(p_1^2\cos \mathrm{\&phi;}-p_2^2\cos \mathrm{\&beta;})-2p_1\sin \mathrm{\&alpha;}\sin \mathrm{\&phi;}\cos \mathrm{\&beta;}+\left[2p_3\cos \mathrm{\&phi;}(\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}-\mathrm{\&rho;}{p}_2)\right]/{\cos }^2\mathrm{\&beta;}}$

In upper two formulas

$\sin \mathrm{\&beta;}=\sin \left(2\mathrm{\&alpha;}\right)\sqrt{1/p^2-{\sin }^2\mathrm{\&phi;}}-\cos \left(2\mathrm{\&alpha;}\right)\sin \mathrm{\&phi;},$ $\cos \mathrm{\&beta;}=\sqrt{1-{\sin }^2\mathrm{\&beta;}},$ $\cos \mathrm{\&phi;}=\sqrt{1-{\sin }^2\mathrm{\&phi;}}$

p ₁＝cosαcosφ+sinαsinφ/ρ，p ₂=cosαcosβ+sinαsinβ/ρ

$p_3=p_1\cos \mathrm{\&phi;}[\frac{\mathrm{\&rho;}\sin \left(2\mathrm{\&alpha;}\right)\sin \mathrm{\&phi;}\cos \mathrm{\&phi;}}{\sqrt{1-{\mathrm{\&rho;}}^2{\sin }^2}\mathrm{\&phi;}}+\cos \left(2\mathrm{\&alpha;}\right)\cos \mathrm{\&phi;}]$

The equal dimensionless of each variable in formula, the span of three parameters is defined as-0.9≤sin φ≤0.9,1,0≤ρ≤2.0 and-0.7≤sin α≤0.7;

To shot point x _swith imaging point (x ₀, T), make V _rmsfor the migration velocity at imaging point place, calculate

$\sin \mathrm{\&phi;}=\frac{x_0-x_s}{\sqrt{{(x_0-x_s)}^2+T^2{V}_{\mathrm{rms}}^2}}$

X in formula ₀and x _sunit be rice, when T is outward journey, unit second, V _rmsunit be meter per second.By rounding sin φ, imaging point place ρ and sin α to determining deviation of three the dimension directions shown, in table, obtain a according to the integer obtained ₁and a ₂; Make f _dbe the dominant frequency of seismic data, unit hertz, calculates

$c^{\&PlusMinus;}=-a_1\&PlusMinus;\sqrt{a_1^2-a_2/\left(2T{f}_d\right)}$

And then calculate two coordinate figures:

$x_r^{\&PlusMinus;}=x+c^{\&PlusMinus;}{\mathrm{TV}}_{\mathrm{rms}}\cos \mathrm{\&alpha;}-{\mathrm{TV}}_{\mathrm{rms}}\frac{\cos \mathrm{\&beta;}(\sin \mathrm{\&beta;}-p_3{c}^{\&PlusMinus;})(1-c^{\&PlusMinus;}\sin \mathrm{\&alpha;}/\mathrm{\&rho;})}{{\cos }^2\mathrm{\&beta;}+p_3\sin \mathrm{\&beta;}{c}^{\&PlusMinus;}-0.5p_3^2{\left(c^{\&PlusMinus;}\right)}^2}$

C in formula ^±dimensionless, unit be rice; Be calculated as follows shot point x again _s, imaging point (x ₀, T-1/ (2f _d)) and acceptance point corresponding incidence and reflection wave angle:

$\cos 2{\mathrm{\&phi;}}_{\mathrm{1,2}}=\frac{1}{\sqrt{1+p_s}\sqrt{1+p_r^{\&PlusMinus;}}}(1+0.5(p_s+p_r^{\&PlusMinus;}-\frac{{(x_s-x_r^{\&PlusMinus;})}^2}{{(T-1/\left(2f_d\right))}^2{V}_{\mathrm{rms}}^2}))$

In formula

$p_s=\frac{{(x_0-x_s)}^2}{V_{\mathrm{rms}}^2{T}^2},$ $p_r^{\&PlusMinus;}=\frac{{(x_0-x_r^{\&PlusMinus;})}^2}{V_{\mathrm{rms}}^2{T}^2}$

Just obtain imaging point (x ₀, T) and corresponding angle domain Fresnel zone [φ ₁, φ ₂], its unit degree of being.

(9) balancing energy is done to the two-dimensional migration road collection of each CDP point, first determine different time degree of depth upper angle territory Fresnel zone, then on angular deflection road collection, excise Fresnel zone part, the remainder on Fresnel zone both sides is superposed respectively, obtain two imaging roads at each CDP point.Fig. 5 is the result of Fig. 4 angular deflection road collection excision Fresnel zone.Fig. 6 is the comparison of stack result before and after Fig. 4 angular deflection road collection excision Fresnel zone.

Automatic gain is done to two-dimensional migration Dao Jizhu road, realizes balancing energy; If the horizontal coordinate of CDP point is x ₀, the shot point coordinate of road centralized shotpoint dimension correspondence is x _s, time depth is T, α is the angle this CDP point different time degree of depth offseting lineups and surface level, take counter clockwise direction as positive-angle, V _rms(T) be the migration velocity in this CDP point different time degree of depth, Δ T is the sampling interval of time depth, calculates

$\sin \mathrm{\&phi;}=\frac{x_0-x_s}{\sqrt{{(x_0-x_s)}^2+T^2{V}_{\mathrm{rms}}^2\left(T\right)}},$ $\mathrm{\&rho;}=\frac{1}{\sqrt{\mathrm{\&Delta;T}}}\sqrt{T-(T-\mathrm{\&Delta;T})\frac{V_{\mathrm{rms}}^2(T-\mathrm{\&Delta;T})}{V_{\mathrm{rms}}^2\left(T\right)}}$

By setting up three-dimensional table spacing used in step 8, sin φ, ρ and sin α being rounded, in table, picking up two real numbers according to the integer obtained, and then calculate the angle domain Fresnel zone [φ in the different time degree of depth by the method for step 8 ₁, φ ₂], its unit degree of being;

In the different time degree of depth, the weight coefficient defined with angle change is

W(θ)＝1.0,θ≤φ ₁-5

$W\left(\mathrm{\&theta;}\right)=0.2+0.8\cos \left(\frac{\mathrm{\&pi;}(\mathrm{\&theta;}-{\mathrm{\&phi;}}_1)}{10}\right),$ φ ₁-5<θ<φ ₁

W(θ)＝0.0,φ ₁≤θ≤φ ₂

$W\left(\mathrm{\&theta;}\right)=0.2+0.8\cos \left(\frac{\mathrm{\&pi;}(\mathrm{\&theta;}-{\mathrm{\&phi;}}_2)}{10}\right),$ φ ₂<θ<φ ₂+5

W(θ)＝1.0,θ≥φ ₂+5

In above formula, θ is that skew road concentrates angle to tie up corresponding angle, unit degree.The angular-trace gather that each shot point is corresponding is concentrated to two-dimensional migration road, each time depth is multiplied by weight coefficient W (θ), just achieve excision Fresnel zone part.Angle corresponding to whole shot point is concentrated in two-dimensional migration road to be greater than and to be less than (φ ₁+ φ ₁the road of)/2 superposes respectively, just obtains two imaging roads at each CDP point.

(10) on the basis at the lineups inclination angle that step 6 obtains, whether the scan method increased and decreased by number percent, completely eliminated reflection line-ups according to scattering wave migration stack section, determine to offset lineups inclination angle accurately.

On the migration stack section obtain step 6, the angle of lineups and surface level carries out number percent increase and decrease, to the universe angle value applying step 9 that each number percent increases and decreases, and by two imaging trace-stackings that each CDP point place obtains, forms stacked section; The migration stack Profile Correlation this section and step 5 obtained, if disappear at a certain imaging point reflection line-ups, angular values is now exactly these lineups inclination angles accurately.

(11) the accurate skew lineups inclination angle utilizing step 10 to obtain, repeat step 9, obtain two migration stack sections be made up of imaging road on the right side of imaging road on the left of Fresnel zone and Fresnel zone respectively, this two section applied energies equilibrium is subtracted each other, obtains only to the scattering wave migration stack section of the independent imaging of little turn-off breakpoint.Fig. 7 and Fig. 8 is the migration stack section that on the left of Fresnel zone, on the right side of imaging road and Fresnel zone, imaging road is formed respectively.Fig. 9 is only to the scattering wave migration stack section of the independent imaging of little turn-off breakpoint.

Utilize reflecting interface and horizontal plane angle accurately, namely lineups inclination angle is offset, repeat step 9, combined respectively in the imaging road on the left of each CDP point Fresnel zone and right side, two the migration imaging sections obtained, remember that its amplitude is A (x respectively, and B (x T), T), wherein x is the horizontal coordinate of CDP point, and T is time depth; Right | A (x, T) | with | B (x, T) | carry out large scale smoothly, obtain its low frequency component, be designated as E (x, T) and F (x, T); To the scattering wave migration stack section of the independent imaging of little turn-off breakpoint be then only

$I(x,T)=(\frac{A(x,T)}{E(x,T)+\mathrm{\&epsiv;}}-\frac{B(x,T)}{F(x,T)+\mathrm{\&epsiv;}})\sqrt{E(x,T)F(x,T)}$

In above formula, ε is a stable factor a small amount of.

(12) by software for display, scattering wave migration stack section numerical value is converted to the profile image of the little breakpoint distribution of subsurface geological structure, this profile image is the tangent plane picture of structure along line direction, this image can indicate the impalpable turn-off of typical reflection ripple migrated image to be less than minor fault and the breakpoint of 1/4 wavelength, for being familiar with migration and the shutoff situation of underground oil and gas, determine effective oil and gas reservoir.Figure 10, Figure 11 and Figure 12 sets forth the little breakpoint distributed image of Fig. 9 and the local contrast of Fig. 3 conventional migration technique stacked section.From Figure 10 and Figure 11, occur the position of little change at reflection line-ups, the imaging of breakpoint scattering wave indicates the existence of breakpoint, and this shows the certain identifiable design breakpoint of the present invention.In Figure 12, the change of lineups is not obvious, and the imaging of breakpoint scattering wave but indicates the existence of breakpoint, and this shows the impalpable little turn-off breakpoint of this anti-bright identifiable design typical reflection ripple migrated image.The final sumbission of this process is: the method can indicate minor fault, the breakpoint of underground structure subtly, thus can determine effective oil and gas reservoir better.

Claims (6)

1. identify a scattering wave Prestack Imaging method for little turn-off breakpoint, it is characterized in that: comprise the following steps: the seismic reflection signals of the artificial epicenter excitation of acquisition mode record of steps A, use wall scroll towing cable or survey line, is recorded on tape; Step B, read seismic signal from tape, form Prestack seismic data, conventional sound attenuation process is done to Prestack seismic data; Build a two dimensional cross-section by along line direction and depth direction, two dimensional processing methods is used to seismic data; For part common midpoint, extract common midpoint gather, conventional NMO velocity pickup is done to the common midpoint gather extracted, the reciprocal value of obtained speed is done laterally average, using the laterally uniform velocity field that obtains as initial offset speed; Step C, foundation initial offset speed, carry out calculations of offset by conventional prestack time method, form CRP gather; Step D, to CRP gather, utilize initial offset speed to do inverse dynamic correction, then do normal moveout correction and obtain the speed after upgrading, do space smoothing to the reciprocal value of the speed after the renewal of each CDP point, namely the level and smooth velocity field obtained is migration velocity field; Step e, foundation migration velocity field, do conventional pre-stack time migration to Prestack seismic data and calculate, and does the dynamic school of high-order residue and stretching excision to the CRP gather obtained, record residual normal moveout and the start offset distance of excising that stretches; Along offset distance superposition CRP gather, obtain migration stack section; Step F, foundation migration stack section and migration velocity field, determining the angle of lineups and surface level on migration stack section, take counter clockwise direction as positive-angle; Step G, foundation migration velocity field, do big gun territory pre-stack time migration to Prestack seismic data, forms the two-dimensional migration road collection with shot point and angle change at each CDP point; The angle of step H, foundation migration velocity field, shot point coordinate, imaging point coordinate, imaging point place skew lineups and surface level, set up the three-dimensional table determining Fresnel zone, three-dimensional table determines the angle domain Fresnel zone of each imaging point accordingly; Step I, balancing energy is done to the two-dimensional migration road collection of each CDP point, first determine different time degree of depth upper angle territory Fresnel zone, then on angular deflection road collection, excise Fresnel zone part, the remainder on Fresnel zone both sides is superposed respectively, obtain two imaging roads at each CDP point; Step J, on the basis at the lineups inclination angle that step F obtains, whether the scan method increased and decreased by number percent, eliminated reflection line-ups completely according on scattering wave migration stack section, determine to offset lineups inclination angle accurately; Step K, the accurate skew lineups inclination angle utilizing step J to obtain, repeat step I, obtain two migration stack sections be made up of imaging road on the right side of imaging road on the left of Fresnel zone and Fresnel zone respectively, this two section applied energies equilibrium is subtracted each other, obtains only to the scattering wave migration stack section of the independent imaging of little turn-off breakpoint; Step L, by software for display scattering wave migration stack section numerical value is converted to the profile image of the little breakpoint distribution of subsurface geological structure, this profile image is the tangent plane picture of structure along line direction, this image can indicate the impalpable turn-off of typical reflection ripple migrated image to be less than minor fault and the breakpoint of 1/4 wavelength, for being familiar with migration and the shutoff situation of underground oil and gas, determine effective oil and gas reservoir.

2. a kind of scattering wave Prestack Imaging method identifying little turn-off breakpoint according to claim 1, it is characterized in that: in step G, described foundation migration velocity field, big gun territory pre-stack time migration is done to Prestack seismic data, the two-dimensional migration road collection formed with shot point and angle change at each CDP point is achieved in that each CDP point, define a three-dimensional array and deposit migration result, the corresponding time depth of one dimension, another dimension is shot point numbering, the third dimension is angle, and angle is from bearing 60 degree to positive 60 degree of changes; Prestack seismic data is pressed common-shot-gather sequence; To the common-shot-gather of monolateral record, utilize big gun, the cautious principle such as mutual of exchanging, the common-shot-gather of the bilateral record of constructing virtual; By shot point coordinate, common-shot-gather is divided into groups, difference group common-shot-gather is placed on the different computing nodes of cluster computer; At each computing node, whole seismic trace is circulated; To each seismic trace, to the whole CDP dot cycles in migration aperture; To each CDP point, determine initial imaging point by migration aperture, the imaging point below this imaging point is calculated successively; First by shot point, geophone station and imaging point coordinate and migration velocity field, when determining to walk and imaging weight coefficient, and then determined to offset amplitude A by the seismologic record amplitude of this seismic trace; Then according to the offset distance of seismic trace, the horizontal coordinate of imaging point and time depth, the residual normal moveout of pickup step e record and the start offset distance of the excision that stretches; The time depth of foundation offset distance and imaging point, judge whether skew amplitude A should be cut, if do not excise, then the angle calculating shot point incident wave and reflecting interface normal is

${\cos }^2\mathrm{\&theta;}=\frac{0.5}{\sqrt{1+p_s}\sqrt{1+p_r}}(1+0.5(p_s+p_r-\frac{H^2}{T^2{V}_{\mathrm{rms}}^2}))+0.5$

In above formula, θ is the angle of shot point incident wave and reflecting interface normal, unit degree, and T is the time depth of imaging point, namely during outward journey, and unit second, V _rmsfor imaging point (x ₀, T) and the migration velocity at place, unit meter per second, H is the offset distance of seismic trace, unit rice, p _sand p _rfor

$p_s=\frac{{(x_0-x_s)}^2}{V_{\mathrm{rms}}^2{T}^2},p_r=\frac{{(x_0-x_g)}^2}{V_{\mathrm{rms}}^2{T}^2}$

Wherein x ₀, x _gand x _sthe horizontal coordinate of imaging point, geophone station and shot point respectively, unit rice;

If Δ θ is the angular interval of predefined, unit degree,

According to cos ²(k Δ θ-Δ θ/2)≤cos ²θ <cos ²(k Δ θ-Δ θ/2) determine positive integer k; If the big gun belonging to this seismic trace is numbered n, if the residual normal moveout of pickup is τ, its unit is second, works as x _g>=x _s, then skew amplitude A is added on (T+ τ, n, the k Δ θ) element of the three-dimensional array depositing migration result, namely in positive-angle; Work as x _g<x _s, then skew amplitude A is added on (T+ τ, the n ,-k Δ θ) element of the three-dimensional array depositing migration result; Complete to whole seismic trace circulation, just obtain the two-dimensional migration road collection with shot point and angle change at each CDP point, this angle is the angle of shot point incident wave and reflecting interface normal.

3. a kind of scattering wave Prestack Imaging method identifying little turn-off breakpoint according to claim 1, it is characterized in that: in steph, described foundation migration velocity field, shot point coordinate, imaging point coordinate, the angle of imaging point place skew lineups and surface level, set up the three-dimensional table determining Fresnel zone, three-dimensional table is determined that angle domain Fresnel zone that each imaging point is corresponding is achieved in that and is made ρ be the ratio of imaging point place actual layer speed and migration velocity accordingly, α is the inclination angle of imaging point place reflecting interface, namely the angle of lineups and surface level is offset, take counter clockwise direction as positive-angle, be that three dimensions set up three-dimensional table with sin φ, ρ and sin α, deposit two dimensionless factors in table, its value is calculated by following two formulas:

$a_1=\frac{\mathrm{\&rho;}(p_1\sin \mathrm{\&phi;}\cos \mathrm{\&beta;}+p_2\sin \mathrm{\&beta;}\cos \mathrm{\&phi;})-\sin \mathrm{\&alpha;}(\cos \mathrm{\&beta;}+\cos \mathrm{\&phi;})}{\mathrm{\&rho;}\cos \mathrm{\&phi;}\cos \mathrm{\&beta;}(p_1^2\cos \mathrm{\&phi;}-p_2^2\cos \mathrm{\&beta;})-2p_1\sin \mathrm{\&alpha;}\sin \mathrm{\&phi;}\cos \mathrm{\&beta;}+\left[2p_3\cos \mathrm{\&phi;}(\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}-\mathrm{\&rho;}{p}_2)\right]/{\cos }^2\mathrm{\&beta;}}$

$a_2=\frac{\mathrm{\&rho;}\cos \mathrm{\&phi;}\cos \mathrm{\&beta;}(\cos \mathrm{\&phi;}+\cos \mathrm{\&beta;})}{\mathrm{\&rho;}\cos \mathrm{\&phi;}\cos \mathrm{\&beta;}(p_1^2\cos \mathrm{\&phi;}-p_2^2\cos \mathrm{\&beta;})-2p_1\sin \mathrm{\&alpha;}\sin \mathrm{\&phi;}\cos \mathrm{\&beta;}+\left[2p_3\cos \mathrm{\&phi;}(\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}-\mathrm{\&rho;}{p}_2)\right]/{\cos }^2\mathrm{\&beta;}}$

In upper two formulas

$\sin \mathrm{\&beta;}=\sin \left(2\mathrm{\&alpha;}\right)\sqrt{1/{\mathrm{\&rho;}}^2-{\sin }^2\mathrm{\&phi;}}-\cos \left(2\mathrm{\&alpha;}\right)\sin \mathrm{\&phi;},\cos \mathrm{\&beta;}=\sqrt{1-{\sin }^2\mathrm{\&beta;}},\cos \mathrm{\&phi;}=\sqrt{1-{\sin }^2\mathrm{\&phi;}}$

p ₁＝cosαcosφ+sinαsinφ/ρ，p ₂＝cosαcosβ+sinαsinβ/ρ

$p_3=p_1\cos \mathrm{\&phi;}[\frac{\mathrm{\&rho;}\sin \left(2\mathrm{\&alpha;}\right)\sin \mathrm{\&phi;}\cos \mathrm{\&phi;}}{\sqrt{1-{\mathrm{\&rho;}}^2{\sin }^2\mathrm{\&phi;}}}+\cos \left(2\mathrm{\&alpha;}\right)\cos \mathrm{\&phi;}]$

The equal dimensionless of each variable in formula, the span of three parameters is defined as-0.9≤sin φ≤0.9,1.0≤ρ≤2.0 and-0.7≤sin α≤0.7;

To shot point x _swith imaging point (x ₀, T), make V _rmsfor the migration velocity at imaging point place, calculate

$\sin \mathrm{\&phi;}=\frac{x_0-x_s}{\sqrt{{(x_0-x_s)}^2+T^2{V}_{\mathrm{rms}}^2}}$

X in formula ₀and x _sunit be rice, when T is outward journey, unit second, V _rmsunit be meter per second; By rounding sin φ, imaging point place ρ and sin α to determining deviation of three the dimension directions shown, in table, obtain a according to the integer obtained ₁and a ₂; Make f _dbe the dominant frequency of seismic data, unit hertz, calculates

$c^{\&PlusMinus;}=-a_1\&PlusMinus;\sqrt{a_1^2-a_2/\left(2{\mathrm{Tf}}_d\right)}$

And then calculate two coordinate figures:

$x_r^{\&PlusMinus;}=x+c^{\&PlusMinus;}{\mathrm{TV}}_{\mathrm{rms}}\cos \mathrm{\&alpha;}-{\mathrm{TV}}_{\mathrm{rms}}\frac{\cos \mathrm{\&beta;}(\sin \mathrm{\&beta;}-p_3{c}^{\&PlusMinus;})(1-c^{\&PlusMinus;}\sin \mathrm{\&alpha;}/\mathrm{\&rho;})}{{\cos }^2\mathrm{\&beta;}+p_3\sin \mathrm{\&beta;}{c}^{\&PlusMinus;}-0.5p_3^2{\left(c^{\&PlusMinus;}\right)}^2}$

C ± dimensionless in formula, unit be rice; Be calculated as follows shot point x again _s, imaging point (x ₀, T-1/ (2f _d)) and acceptance point corresponding incidence and reflection wave angle:

$\cos 2{\mathrm{\&phi;}}_{\mathrm{1,2}}=\frac{1}{\sqrt{1+p_s}\sqrt{1+p_r^{\&PlusMinus;}}}(1+0.5(p_s+p_r^{\&PlusMinus;}-\frac{{(x_s-x_r^{\&PlusMinus;})}^2}{(T-1/{\left(2f_d)\right)}^2{V}_{\mathrm{rms}}^2}))$

In formula

$p_s=\frac{{(x_0-x_s)}^2}{V_{\mathrm{rms}}^2{T}^2},p_r^{\&PlusMinus;}=\frac{{(x_0-x_r^{\&PlusMinus;})}^2}{V_{\mathrm{rms}}^2{T}^2}$

Just obtain imaging point (x ₀, T) and corresponding angle domain Fresnel zone [φ ₁, φ ₂], its unit degree of being.

4. a kind of scattering wave Prestack Imaging method identifying little turn-off breakpoint according to claim 1, it is characterized in that: in step I, described does balancing energy to the two-dimensional migration road collection of each CDP point, first determine different time degree of depth upper angle territory Fresnel zone, then on angular deflection road collection, Fresnel zone part is excised, the remainder on Fresnel zone both sides is superposed respectively, obtain two imaging roads at each CDP point to be achieved in that and to do automatic gain to two-dimensional migration Dao Jizhu road, realize balancing energy; If the horizontal coordinate of CDP point is x ₀, the shot point coordinate of road centralized shotpoint dimension correspondence is x _s, time depth is T, α is the angle this CDP point different time degree of depth offseting lineups and surface level, take counter clockwise direction as positive-angle, V _rms(T) be the migration velocity in this CDP point different time degree of depth, Δ T is the sampling interval of time depth, calculates

$\sin \mathrm{\&phi;}=\frac{x_0-x_s}{\sqrt{{(x_0-x_s)}^2+T^2{V}_{\mathrm{rms}}^2\left(T\right)}},\mathrm{\&rho;}=\frac{1}{\sqrt{\mathrm{\&Delta;T}}}\sqrt{T-(T-\mathrm{\&Delta;T})\frac{V_{\mathrm{rms}}^2(T-\mathrm{\&Delta;T})}{V_{\mathrm{rms}}\left(T\right)}}$

By setting up three-dimensional table spacing used in step H, sin φ, ρ and sin α being rounded, in table, picking up two real numbers according to the integer obtained, and then calculate the angle domain Fresnel zone [φ in the different time degree of depth by the method for step H ₁, φ ₂], its unit degree of being;

In the different time degree of depth, the weight coefficient defined with angle change is

W(θ)＝1.0，θ≤φ ₁-5

$W\left(\mathrm{\&theta;}\right)=0.2+0.8\cos \left(\frac{\mathrm{\&pi;}(\mathrm{\&theta;}-{\mathrm{\&phi;}}_1)}{10}\right),{\mathrm{\&phi;}}_1-5<\mathrm{\&theta;}<{\mathrm{\&phi;}}_1$

W(θ)＝0.0，φ ₁≤θ≤φ ₂

$W\left(\mathrm{\&theta;}\right)=0.2+0.8\cos \left(\frac{\mathrm{\&pi;}(\mathrm{\&theta;}-{\mathrm{\&phi;}}_2)}{10}\right),{\mathrm{\&phi;}}_2<\mathrm{\&theta;}<{\mathrm{\&phi;}}_2+5$

W(θ)＝1.0，θ≥φ ₂+5

In above formula, θ is that skew road concentrates angle to tie up corresponding angle, unit degree; The angular-trace gather that each shot point is corresponding is concentrated to two-dimensional migration road, each time depth is multiplied by weight coefficient W (θ), just achieve excision Fresnel zone part; Angle corresponding to whole shot point is concentrated in two-dimensional migration road to be greater than and to be less than (φ ₁+ φ ₁the road of)/2 superposes respectively, just obtains two imaging roads at each CDP point.

5. a kind of scattering wave Prestack Imaging method identifying little turn-off breakpoint according to claim 1, it is characterized in that: in step J, described on the basis at the lineups inclination angle that step F obtains, by the scan method that number percent increases and decreases, whether reflection line-ups is eliminated completely according on scattering wave migration stack section, determine to offset accurately lineups inclination angle and be achieved in that the angle of lineups and surface level on the migration stack section that obtain step F carries out number percent increase and decrease, to the universe angle value applying step I that each number percent increases and decreases, and by two imaging trace-stackings that each CDP point place obtains, form stacked section, the migration stack Profile Correlation this section and step e obtained, if disappear at a certain imaging point reflection line-ups, angular values is now exactly these lineups inclination angles accurately.

6. a kind of scattering wave Prestack Imaging method identifying little turn-off breakpoint according to claim 1, it is characterized in that: in step K, the described accurate skew lineups inclination angle utilizing step J to obtain, repeat step I, obtain two migration stack sections be made up of imaging road on the right side of imaging road on the left of Fresnel zone and Fresnel zone respectively, this two section applied energies equilibrium is subtracted each other, obtain only being achieved in that the scattering wave migration stack section of the independent imaging of little turn-off breakpoint and utilize reflecting interface and horizontal plane angle accurately, namely lineups inclination angle is offset, repeat step I, the imaging road on the left of each CDP point Fresnel zone and right side is combined respectively, two the migration imaging sections obtained, remember that its amplitude is A (x respectively, and B (x T), T), wherein x is the horizontal coordinate of CDP point, T is time depth, right | A (x, T) | with | B (x, T) | carry out large scale smoothly, obtain its low frequency component, be designated as E (x, T) and F (x, T), to the scattering wave migration stack section of the independent imaging of little turn-off breakpoint be then only

$I(x,T)=(\frac{A(x,T)}{E(x,T)+\mathrm{\&epsiv;}}-\frac{B(x,T)}{F(x,T)+\mathrm{\&epsiv;}})\sqrt{E(x,T)F(x,T)}$

In above formula, ε is a stable factor a small amount of.