CN102866421A - Scattered wave pre-stack imaging method for identifying small-fault throw breakpoints - Google Patents

Abstract

The invention relates to a scattered wave pre-stack imaging method for identifying small-fault throw breakpoints, which is applied to processing of seismic reflection information in seismic exploration and is of a seismic imaging method of a micro-small fault for identifying the underground geological structure. According to the method, scattered waves of the breakpoints of the fault are utilized for independent imaging of the breakpoints with the fault throw of less than 1/4 wavelength, which are difficult to identify by conventional reflected wave migration. According to the method, the existing seismic reflection information is utilized, and the independent imaging of the small breakpoints is realized by pressing the energy of transmitted waves in pre-stack migration calculation and utilizing the polarity inversion characteristic of the scattered waves of the breakpoints. Pre-stack time migration is taken as the basis, an angle gather is generated by migration of single-shot information, and a Fresnel zone corresponding to an inphase axis of reflected waves is excised so as to eliminate the transmitted waves. According to the method, the small fault and breakpoints of the underground structure can be indicated in a fine manner, migration and plugging of underground oil gas can be better known, and the method further has an important application value for exploration and development of oil gas resources. The ability to identify the underground deep fissures and fault of the seismic exploration method is significantly improved.

Description

Identify the scattering wave prestack formation method of little turn-off breakpoint

Technical field

The invention belongs to reflected seismic information processing technology field in the seismic prospecting, relating to the prestack migration image technology category in the process of seismic data processing, is a kind of offset imaging method that can identify the small tomography of tectonic structure that the 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 that understanding oil gas is dredged system, and then the important step of identification Favorable Reservoir.Therefore, in the reflected seismic information processing procedure, the resolution that improves migration imaging is an important effort target always in seismic prospecting.

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

For improving the resolution of seismic migration imaging, develop many methods, comprised the spectral whitening deconvolution, unstable state deconvolution of imaging stacked section, based on all kinds of band technology, anti-Q filtering, the viscoelasticity pre-stack depth migrations etc. widened of statistical hypothesis or well-log information.The core of each class methods all is by improving the upper limit of seismic event frequency band, so that the reflection line-ups of playback becomes is thinner, thereby improves 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 satisfy oil-gas exploration and exploitation to the needs of minor fault identification.

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

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

A kind of method of suppressing the interface echo imaging is to reject interface echo in seismic data, and this is to utilize reflection wave some feature in Prestack seismic data, realizes such as the Hyperbolic Feature in common midpoint gather.But what the feature of interface echo became in the complex geological structure situation is more complicated, and particularly in the lower situation of seismic data signal to noise ratio (S/N ratio), it almost is a job that is difficult to finish for accurate perception reflex ripple lineups and correct rejecting; In addition, adopt the transform method such as radon transform also will destroy the very weak scattering wave of script.Therefore the present invention adopts in the migration imaging road and concentrates the depressor reflex ripple, realizes the scattering wave imaging.

Because scattering wave amplitude and all kinds of noise be in same magnitude, even also hang down a magnitude than some relevant noises.For avoiding noise to the interference of scattering wave imaging results, the present invention utilizes the scattering wave of little breakpoint to have the characteristics of reversal of poles, and namely the polarity of scattering wave lineups is strengthened the imaging to little breakpoint in the anti-living positive and negative conversion of certain point.The present invention is positioned the scattering wave of little turn-off breakpoint is carried out imaging, and these small breakpoints are vital to the identification of hydrocarbon system of dredging.

Pre-stack time migration is a kind of important method in the seismic migration imaging, compares with prestack depth migration method, and except having higher counting yield, its main advantage is only to need to use stack (root mean square) speed; The mode such as Negotiation speed scanning obtains appropriate rate pattern so simply.Because it is complicated that the actual in large quantities Geological Problems that runs into all is reflective construct, but the horizontal change of medium velocity is not very violent, so the present invention is take prestack time migration method as base growth scattering wave prestack migration image method.This has avoided the main difficulty of a migration imaging: velocity modeling.Imaging is even more important this point to scattering wave, because scattering wave is more responsive to the slight error of speed, needs to use the more accurate rate pattern of more conventional pre-stack depth migration.

Summary of the invention

The objective of the invention is: a kind of scattering wave prestack formation method of identifying little turn-off breakpoint is provided, and it can suppress the transmitted wave energy in the migration process of seismic data, obtain the independent imaging to little turn-off breakpoint; This method can improve the reflected seismic information migration imaging to the recognition capability of small tomography, improves the understanding of in the oil-gas exploration and development Deep Oil And Gas Exploration being dredged system.

The technical solution used in the present invention is: identify the scattering wave prestack formation method of little turn-off breakpoint, concrete steps comprise:

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

(2) read seismic signal from tape, form Prestack seismic data, Prestack seismic data is done conventional sound attenuation process; By making up a two dimensional cross-section along line direction and depth direction, seismic data is used the two-dimensional process method; For the part common midpoint, extract common midpoint gather, the common midpoint gather that extracts is made conventional NMO velocity picks up, the reciprocal value of resulting speed is done laterally average, with obtain laterally uniformly velocity field as initial offset speed;

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

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

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

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

(7) according to the migration velocity field, Prestack seismic data is done big gun territory pre-stack time migration, form the two-dimensional migration road collection that changes with shot point and angle 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 of determining Fresnel zone, three-dimensional table is determined the angle domain Fresnel zone of each imaging point accordingly;

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

(10) on the basis at the lineups inclination angle that step 6 obtains, by the scan method of number percent increase and decrease, whether eliminated reflection line-ups fully according to scattering wave skew stacked section, determine to be offset accurately the lineups inclination angle;

(11) the accurate skew lineups inclination angle that utilizes step 10 to obtain, repeating step 9, obtain two skew stacked sections that consisted of by Fresnel zone imaging road, left side and imaging road, Fresnel zone right side respectively, this two section applied energies equilibrium is subtracted each other, only obtain the scattering wave skew stacked section to the independent imaging of little turn-off breakpoint;

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

Described according to the migration velocity field, Prestack seismic data is done big gun territory pre-stack time migration, forming the two-dimensional migration road collection that changes with shot point and angle 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 the shot point numbering, and the third dimension is angle, and angle changes to positive 60 degree from negative 60 degree; Prestack seismic data is pressed the common-shot-gather ordering; 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 the shot point coordinate common-shot-gather is divided into groups, with not on the same group common-shot-gather be placed on the different computing nodes of cluster computer; At each computing node, to whole seismic trace circulations; To each seismic trace, to the whole CDP dot cycles in the migration aperture; To each CDP point, determine initial imaging point by migration aperture, this imaging point below imaging point is calculated successively; First by shot point, geophone station and imaging point coordinate and migration velocity field, when determining to walk and the imaging weight coefficient, and then determine skew amplitude A by the seismologic record amplitude of this seismic trace; Then according to the offset distance of seismic trace, horizontal coordinate and the time depth of imaging point, pick up the start offset distance of the residual normal moveout of step 5 record and the excision that stretches; According to the time depth of offset distance and imaging point, judge whether skew amplitude A should be cut, if do not excise, the angle that then calculates shot point incident wave and reflecting interface normal is

${\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">cos</figure-callout>}}^2\mathrm{\&theta;}=\frac{0.5}{\sqrt{1+{\mathrm{<figure-callout\; id="1"\; label="p\; s"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">p</figure-callout>}}_s}\sqrt{1+p_r}}(1+0.5(p_s+p_r-\frac{H^2}{T^2{V}_{\mathrm{rms}}^2}))+0.5$

In the following formula, θ is the angle of shot point incident wave and reflecting interface normal, unit degree, and T is the time depth of imaging point, i.e. during outward journey, unit second, V _RmsBe imaging point (x ₀, the migration velocity of T) locating, unit meter per second, H are the offset distances of seismic trace, unit rice, p _sAnd p _rFor

$p_s=\frac{{(x_0-x_s)}^2}{V_{\mathrm{rms}}^2{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">T</figure-callout>}}^2},$ $p_r=\frac{{(x_0-x_g)}^2}{V_{\mathrm{rms}}^2{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">T</figure-callout>}}^2}$

X wherein ₀, x _gAnd x _sRespectively the horizontal coordinate of imaging point, geophone station and shot point, unit rice;

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

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 of determining Fresnel zone, three-dimensional table is determined that angle domain Fresnel zone corresponding to each imaging point is achieved in that and is made that ρ is the ratio of imaging point place actual layer speed and migration velocity accordingly, α is the inclination angle of imaging point place reflecting interface, namely be offset the angle of lineups and surface level, take counterclockwise as positive-angle; Set up three-dimensional table take sin φ, ρ and sin α as three dimensions, deposit two dimensionless factors in the table, its value is calculated by following two formulas:

$a_1=\frac{\mathrm{\&rho;}({\mathrm{<figure-callout\; id="1"\; label="p"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">p</figure-callout>}}_1\sin \mathrm{\&phi;}\cos \mathrm{\&beta;}+{\mathrm{<figure-callout\; id="2"\; label="p"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">p</figure-callout>}}_2\sin \mathrm{\&beta;}\cos \mathrm{\&phi;})-\sin \mathrm{\&alpha;}(\cos \mathrm{\&beta;}+\cos \mathrm{\&phi;})}{\mathrm{\&rho;}\cos \mathrm{\&phi;}\cos \mathrm{\&beta;}({\mathrm{<figure-callout\; id="1"\; label="p"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">p</figure-callout>}}_1^2\cos \mathrm{\&phi;}-{\mathrm{<figure-callout\; id="2"\; label="p"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">p</figure-callout>}}_2^2\cos \mathrm{\&beta;})-2{\mathrm{<figure-callout\; id="1"\; label="p"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">p</figure-callout>}}_1\sin \mathrm{\&alpha;}\sin \mathrm{\&phi;}\cos \mathrm{\&beta;}+\left[2{\mathrm{<figure-callout\; id="3"\; label="p"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">p</figure-callout>}}_3\cos \mathrm{\&phi;}(\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}-\mathrm{\&rho;}{p}_2)\right]/{\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">cos</figure-callout>}}^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;}({\mathrm{<figure-callout\; id="1"\; label="p"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">p</figure-callout>}}_1^2\cos \mathrm{\&phi;}-{\mathrm{<figure-callout\; id="2"\; label="p"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">p</figure-callout>}}_2^2\cos \mathrm{\&beta;})-2{\mathrm{<figure-callout\; id="1"\; label="p"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">p</figure-callout>}}_1\sin \mathrm{\&alpha;}\sin \mathrm{\&phi;}\cos \mathrm{\&beta;}+\left[2{\mathrm{<figure-callout\; id="3"\; label="p"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">p</figure-callout>}}_3\cos \mathrm{\&phi;}(\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}-\mathrm{\&rho;}{p}_2)\right]/{\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">cos</figure-callout>}}^2\mathrm{\&beta;}}$

In upper two formulas

$\sin \mathrm{\&beta;}=\sin \left(2\mathrm{\&alpha;}\right)\sqrt{1/p^2-{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">sin</figure-callout>}}^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-{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2\mathrm{\&phi;}}$

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

${\mathrm{<figure-callout\; id="3"\; label="p"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">p</figure-callout>}}_3={\mathrm{<figure-callout\; id="1"\; label="p"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">p</figure-callout>}}_1\cos \mathrm{\&phi;}[\frac{\mathrm{\&rho;}\sin \left(2\mathrm{\&alpha;}\right)\sin \mathrm{\&phi;}\cos \mathrm{\&phi;}}{\sqrt{1-{\mathrm{\&rho;}}^2{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2}\mathrm{\&phi;}}+\cos \left(2\mathrm{\&alpha;}\right)\cos \mathrm{\&phi;}]$

The equal dimensionless of each variable in the 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 _sAnd imaging point (x ₀, T), make V _RmsBe the migration velocity at imaging point place, calculate

$\sin \mathrm{\&phi;}=\frac{x_0-x_s}{\sqrt{{(x_0-x_s)}^2+T^2{\mathrm{<figure-callout\; id="2"\; label="V\; rms"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">V</figure-callout>}}_{\mathrm{rms}}^{}}}$

X in the formula ₀And x _sUnit be rice, when T is outward journey, unit second, V _RmsUnit be meter per second.By sin φ, the ρ of imaging point place and sin α being rounded to determining deviation of three dimension directions showing, in table, obtain a according to the integer that obtains ₁And a ₂Make f _dBe the dominant frequency of seismic data, the unit hertz calculates

$c^{\&PlusMinus;}=-{\mathrm{<figure-callout\; id="1"\; label="a"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">a</figure-callout>}}_1\&PlusMinus;\sqrt{{\mathrm{<figure-callout\; id="1"\; label="a"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">a</figure-callout>}}_1^2-{\mathrm{<figure-callout\; id="2"\; label="a"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">a</figure-callout>}}_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;})}{{\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">cos</figure-callout>}}^2\mathrm{\&beta;}+{\mathrm{<figure-callout\; id="3"\; label="p"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">p</figure-callout>}}_3\sin \mathrm{\&beta;}{c}^{\&PlusMinus;}-0.5{\mathrm{<figure-callout\; id="3"\; label="p"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">p</figure-callout>}}_3^2{\left(c^{\&PlusMinus;}\right)}^2}$

C in the formula ^±Dimensionless,

Unit be rice; Be calculated as follows again shot point x _s, imaging point (x ₀, T-1/ (2f _d)) and acceptance point

Corresponding incident and reflection wave angle:

$\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">cos</figure-callout>}2{\mathrm{\&phi;}}_{\mathrm{1,2}}=\frac{1}{\sqrt{1+{\mathrm{<figure-callout\; id="1"\; label="p\; s"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">p</figure-callout>}}_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 the formula

$p_s=\frac{{(x_0-x_s)}^2}{V_{\mathrm{rms}}^2{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">T</figure-callout>}}^2},$ $p_r^{\&PlusMinus;}=\frac{{(x_0-x_r^{\&PlusMinus;})}^2}{V_{\mathrm{rms}}^2{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">T</figure-callout>}}^2}$

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

The described two-dimensional migration road collection that each CDP is ordered is done balancing energy, determine first different time degree of depth upper angle territory Fresnel zone, then in angular deflection road collection excision Fresnel zone part, the remainder on Fresnel zone both sides is superposeed respectively, obtain two imaging roads at each CDP point and be achieved in that two-dimensional migration road collection is pursued the road does automatic gain, realizes balancing energy; If the horizontal coordinate that CDP is ordered is x ₀, shot point coordinate corresponding to road centralized shotpoint dimension is x _s, time depth is T, α is for the angle of skew lineups and surface level on this CDP point different time degree of depth, and take counter clockwise direction as positive-angle, V _Rms(T) be the migration velocity on 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 the used spacing of three-dimensional table in the step 8 sin φ, ρ and sin α are rounded, in table, pick up two real numbers according to the integer that obtains, and then by the angle domain Fresnel zone [φ on the method calculating different time degree of depth of step 8 ₁, φ ₂], its unit degree of being;

In the different time degree of depth, definition with the weight coefficient that angle changes 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

θ concentrates angle corresponding to angle dimension, unit degree for the skew road in the following formula.Angular-trace gather corresponding to each shot point concentrated in the two-dimensional migration road, be multiply by weight coefficient W (θ) at each time depth, just realized excision Fresnel zone part.With the two-dimensional migration road concentrate angle corresponding to whole shot points greater than with 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, scan method by the number percent increase and decrease, according to whether having rejected reflection line-ups fully on the scattering wave skew stacked section, determine to be offset accurately the lineups inclination angle and be achieved in that on the skew stacked section that step 6 is obtained that the angle of lineups and surface level carries out number percent and increases and decreases, universe angle value applying step 9 to each number percent increase and decrease, and with two imaging trace-stackings that each CDP point place obtains, form stacked section; With the skew stacked section contrast that this section and step 5 obtain, if disappear at a certain imaging point reflection line-ups, the angle numerical value of this moment is exactly accurately inclination angles of these lineups.

The described accurate skew lineups inclination angle that utilizes step 10 to obtain, repeating step 9, obtain two skew stacked sections that consisted of by Fresnel zone imaging road, left side and imaging road, Fresnel zone right side respectively, this two section applied energies equilibrium is subtracted each other, only obtaining scattering wave skew stacked section to the independent imaging of little turn-off breakpoint is achieved in that and utilizes accurately reflecting interface and horizontal plane angle, namely be offset the lineups inclination angle, repeating step 9, the imaging road on each CDP point Fresnel zone left side and right side is made up respectively, two migration imaging sections that obtain, remember that respectively its amplitude is A (x, T) and B (x, T), wherein x is the horizontal coordinate that CDP is ordered, and T is time depth; Right | A (x, T) | and | B (x, T) | it is level and smooth to carry out large scale, obtains its low frequency component, is designated as E (x, T) and F (x, T); Only to the scattering wave skew stacked section of the independent imaging of little turn-off breakpoint then be

$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)}$

ε is a stable factor a small amount of in the following formula.

The scattering wave prestack formation 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 formation method of the little turn-off breakpoint of identification of the present invention can to the impalpable turn-off of conventional reflection wave imaging less than the independent imaging of the breakpoint of 1/4 wavelength, can significantly improve method of seismic prospecting to the recognition capability of underground crack and tomography.

The scattering wave prestack formation method of the little turn-off breakpoint of identification of the present invention can be suppressed the transmitted wave energy in migration process, outstanding breakpoint scattering wave.

The scattering wave prestack formation method of the little turn-off breakpoint of identification of the present invention can independently determine migration velocity and interface dip that the scattering wave skew is used.

Specific implementation principle of the present invention is as follows:

The present invention adopts in the migration imaging road and concentrates the depressor reflex ripple to realize the scattering wave imaging, and utilizes the reversal of poles of little breakpoint scattering wave, strengthens little breakpoint imaging; Its core is the different manifestations on single big gun skew road collection according to reflection line-ups and scattering wave lineups, comes the depressor reflex ripple, keeps scattering wave.Mainly contain following 3 points, the one, improve prestack time migration method, obtain the angle domain skew road collection of big gun record; The 2nd, the accurate Fresnel zone that determines reflection wave on angle domain skew road collection is taken this depressor reflex ripple; The 3rd, breakpoint scattering wave image enhancement.

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

1. angle domain is offset the road collection

The prestack time can be regarded pre-stack depth migration as in the velocity equivalent parameter, i.e. being similar under the root-mean-square velocity.Phase-shift method and steady phase point principle based on pre-stack depth migration, can derive conventional pre-stack time migration formula, and based on this framework, can be in the situation that does not need the known interface normal, obtain the angle of incident wave and interface normal, take this to obtain the angle domain skew road collection of incident wave corresponding to single shot record and interface normal angle.

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 regard single seismic trace as a big gun record, further nonhomogeneous media is approximately layered medium, by the phase-shift method of the degree of depth-wave number-frequency field 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{{\mathrm{<figure-callout\; id="2"\; label="v\; i"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">v</figure-callout>}}_i^{}}{{\mathrm{\&omega;}}^2}{k}_x^2}\right)---\left(1\right)$

Express the degree of depth with the straight whilst on tour T of one way that hangs down in the formula,

Δ T _iBe the thickness that each layer medium expressed during with outward journey, Δ z is arranged _i=v _iΔ T _i, v _iBe the speed of each layer medium, n is the above medium number of plies of zone of interest,

It is the wave number-frequency field continuation result at time depth T place.

Formula (1) but cumulative part approximate expression in the exponential term of Exponential function be:

$\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left(\mathrm{\&Delta;}{T}_i\sqrt{1-k_x^2{\mathrm{<figure-callout\; id="2"\; label="v\; i"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">v</figure-callout>}}_i^{}/{\mathrm{\&omega;}}^2}\right)\&ap;\sqrt{1-k_x^2{\mathrm{<figure-callout\; id="2"\; label="V\; rms"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">V</figure-callout>}}_{\mathrm{rms}}^{}/{\mathrm{\&omega;}}^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:

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

V _RmsNamely be the needed migration velocity of pre-stack time migration, i.e. conventional root-mean-square velocity.

With formula (2) substitution (1) formula and do the space inverse fourier transform, get the space of single seismic trace-frequency field continuation result and be:

$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)$

Introduced ray parameter p in the formula _x=k _x/ ω.The available steady phase point principle of formula (4) is tried to achieve progressive solution, and steady ray parameter corresponding to phase point then represented 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)$

With the steady phase point principle of formula (6) substitution formula, 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;}}_{\mathrm{<figure-callout\; id="3"\; label="g"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">g</figure-callout>}}^{3/2}{V}_{\mathrm{rms}}}\exp \left(\mathrm{j\&omega;}{\mathrm{\&tau;}}_g\right)---\left(7\right)$

τ in the formula _gBe seimic travel time, have

${\mathrm{\&tau;}}_g=\mathrm{<figure-callout\; id="1"\; label="T"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">T</figure-callout>}\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 be considered the propagation of seismic event in three dimensions, so the amplitude item in the 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{\mathrm{<figure-callout\; id="2"\; label="V\; rms"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">V</figure-callout>}}_{\mathrm{rms}}^{}}\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 in the different time degree of depth is

$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{\mathrm{<figure-callout\; id="2"\; label="V\; rms"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">V</figure-callout>}}_{\mathrm{rms}}^{}}\exp (-j{\mathrm{\&omega;\&tau;}}_s)---\left(11\right)$

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

${\mathrm{\&tau;}}_s=\mathrm{<figure-callout\; id="1"\; label="T"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">T</figure-callout>}\sqrt{1+\frac{{(x-x_s)}^2}{V_{\mathrm{rms}}^2{T}^2}}---\left(12\right)$

Source wavelet is unknown generally speaking.Since deconvolution can be rejected the impact of wavelet, we can ignore the impact of source wavelet continuous item in imaging, namely

Record n seismic trace if establish big gun, use the deconvolution image-forming condition, the imaging results that can get 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 the formula It is half derivative of seismic trace time-domain signal; Be the acceptance point of m seismic trace during to the walking of imaging point, can obtain by the horizontal coordinate substitution formula (8) with this acceptance point.

To locate the interval velocity of medium be v if assume picture point (x, T), and according to the ray parameter of formula (6) and (10), the remaining device for carrying a tripot of wavefront surface direction that can obtain shot point incident field and acceptance point anti-pass wave field is vp _xWith

And then can obtain shot point incident wave corresponding to this seismic trace and acceptance point anti-pass ripple at the angle of imaging point is

$\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">cos</figure-callout>}2\mathrm{\&theta;}=\frac{1}{\sqrt{1+p_s}\sqrt{1+p_r}}(\sqrt{1+(1-{\mathrm{\&rho;}}^2){\mathrm{<figure-callout\; id="1"\; label="p\; s"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">p</figure-callout>}}_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)$

2 θ are angles of incident wave and reflection wave in the formula, ρ=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{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">T</figure-callout>}}^2},$ $p_r=\frac{{(x-x_g)}^2}{V_{\mathrm{rms}}^2{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">T</figure-callout>}}^2}$

For distinguishing better reflection and scattering wave, the angle that the present invention tries to achieve formula (14) is introduced positive and negative, definition x _g〉=x _sThe time angle for just.

Practical application middle level speed v is unknown.Because angle Dao Ji of the present invention only is used for the depressor reflex ripple, rather than utilizes angle-magnitude relation to carry out inverting, therefore when definite angle θ, can suppose ρ=1, this pattern (14) but approximate expression be

$\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">cos</figure-callout>}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 the angle varied pitch, in formula (13) cumulative to each imaging point (x, T) consider the angle that each seismic trace is corresponding, it is cumulative respectively to each angular range to press k Δ θ-Δ θ/2≤θ＜k Δ θ+Δ θ/2, just can obtain the angle domain skew road collection I (x of single shot record, T, k Δ θ).For reducing the calculated amount of directly being calculated angle by formula (15), can obtain cos by formula (15) ²θ is with itself and cos ²(k Δ θ ± Δ θ/2) relatively determines the cumulative position of this skew amplitude.Work as x _r〉=x _s, the skew amplitude is added on the positive-angle, otherwise, then be added to negative angle.

2. angle domain Fresnel zone

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

So-called Fresnel zone namely is in this zone, and the time depth of the skew lineups of reflection wave and imaging point place skew lineups is poor less than specific spacing, for example 1/4 cycle; Because seismic signal is in a frequency band, the present invention is defined as 1/2 cycle corresponding to seismic event dominant frequency with this time interval.The below will based on the ray parameter expression formula of formula (6) and (10), try to achieve the angle of Fresnel zone two ends correspondence on Jiao Daoji at inclined reflection interface, i.e. the angle domain Fresnel zone.

As shown in Figure 1, the inclination angle of establishing reflecting interface is α, and the coordinate of shot point S is x _s, the coordinate of imaging point I is (x, T), it also is the incidence point of incident wave on reflecting interface; If on the interface, be that the point of d is the border of Fresnel zone along the interface apart from imaging point, i.e. I ' point, then the picture concentrated in road, imaging point angle of this point reflection ripple will be located at (x, T-δ), 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 and the angle identical (intersect at I ' point) of reflection wave with interface normal, and shot point arrives x again to Fresnel zone frontier point (I ' point) _rThe point walk the time with shot point to the P point again to R order walk the time identical these two conditions, 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 the formula _Rms, d is the distance along direction imaging point Fresnel zone border, interface, x _rIt is the acceptance point coordinate that plan is found the solution.Solve two groups of (x by formula (16) and (17) _r, d), with x _s, x _r(x, T-1/ (2f _d)) substitution formula (15) can try to achieve Fresnel zone two ends corresponding angle on Jiao Daoji.

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

${\mathrm{<figure-callout\; id="2"\; label="v"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">v</figure-callout>}}^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})]$

Thereby

${\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)$

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

This group nonlinear equation of direct solution formula (16) and (17) is difficult, utilizes d and TV in the formula (16) _RmsComparing is in a small amount, and 1/ (2f in the formula (17) _d) to compare with T be in a small amount, can be 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 take sin φ, ρ and sin α as dimension, deposit two real numbers of trying to achieve by as shown in the formula (19) and (20) in the table:

$a_1=\frac{\mathrm{\&rho;}({\mathrm{<figure-callout\; id="1"\; label="p"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">p</figure-callout>}}_1\sin \mathrm{\&phi;}\cos \mathrm{\&beta;}+{\mathrm{<figure-callout\; id="2"\; label="p"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">p</figure-callout>}}_2\sin \mathrm{\&beta;}\cos \mathrm{\&phi;})-\sin \mathrm{\&alpha;}(\cos \mathrm{\&beta;}+\cos \mathrm{\&phi;})}{\mathrm{\&rho;}\cos \mathrm{\&phi;}\cos \mathrm{\&beta;}({\mathrm{<figure-callout\; id="1"\; label="p"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">p</figure-callout>}}_1^2\cos \mathrm{\&phi;}-{\mathrm{<figure-callout\; id="2"\; label="p"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">p</figure-callout>}}_2^2\cos \mathrm{\&beta;})-2{\mathrm{<figure-callout\; id="1"\; label="p"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">p</figure-callout>}}_1\sin \mathrm{\&alpha;}\sin \mathrm{\&phi;}\cos \mathrm{\&beta;}+\left[2{\mathrm{<figure-callout\; id="3"\; label="p"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">p</figure-callout>}}_3\cos \mathrm{\&phi;}(\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}-\mathrm{\&rho;}{p}_2)\right]/{\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">cos</figure-callout>}}^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;}({\mathrm{<figure-callout\; id="1"\; label="p"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">p</figure-callout>}}_1^2\cos \mathrm{\&phi;}-{\mathrm{<figure-callout\; id="2"\; label="p"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">p</figure-callout>}}_2^2\cos \mathrm{\&beta;})-2{\mathrm{<figure-callout\; id="1"\; label="p"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">p</figure-callout>}}_1\sin \mathrm{\&alpha;}\sin \mathrm{\&phi;}\cos \mathrm{\&beta;}+\left[2{\mathrm{<figure-callout\; id="3"\; label="p"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">p</figure-callout>}}_3\cos \mathrm{\&phi;}(\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}-\mathrm{\&rho;}{p}_2)\right]/{\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">cos</figure-callout>}}^2\mathrm{\&beta;}}---\left(20\right)$

In upper two formulas

$\sin \mathrm{\&beta;}=\sin \left(2\mathrm{\&alpha;}\right)\sqrt{1/p^2-{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2\mathrm{\&phi;}}-\cos \left(2\mathrm{\&alpha;}\right)\sin \mathrm{\&phi;},$ $\cos \mathrm{\&beta;}=\sqrt{1-{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2\mathrm{\&beta;}},$ $\cos \mathrm{\&phi;}=\sqrt{1-{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2\mathrm{\&phi;}}$

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

${\mathrm{<figure-callout\; id="3"\; label="p"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">p</figure-callout>}}_3={\mathrm{<figure-callout\; id="1"\; label="p"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">p</figure-callout>}}_1\cos \mathrm{\&phi;}[\frac{\mathrm{\&rho;}\sin \left(2\mathrm{\&alpha;}\right)\sin \mathrm{\&phi;}\cos \mathrm{\&phi;}}{\sqrt{1-{\mathrm{\&rho;}}^2{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">sin</figure-callout>}}^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{\mathrm{<figure-callout\; id="2"\; label="V\; rms"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">V</figure-callout>}}_{\mathrm{rms}}^{}}}$

ρ and sin α corresponding to inclination alpha according to sin φ, imaging point place pick up a in table ₁And a ₂Calculate

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

Can obtain x _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;})}{{\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">cos</figure-callout>}}^2\mathrm{\&beta;}+{\mathrm{<figure-callout\; id="3"\; label="p"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">p</figure-callout>}}_3\sin \mathrm{\&beta;}{c}^{\&PlusMinus;}-0.5{\mathrm{<figure-callout\; id="3"\; label="p"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">p</figure-callout>}}_3^2{\left(c^{\&PlusMinus;}\right)}^2}---\left(21\right)$

With formula (21) and imaging point (x, T-1/ (2f _d)) substitution formula (15), can solve two angle φ corresponding to Fresnel zone ₁And φ ₂

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

When determining Fresnel zone, need utilize the inclination angle of reflecting interface, i.e. the angle of reflection line-ups and surface level; Because what prestack time migration method obtained is the time domain migrated section, generally only can obtain approximate reflector dip based on this section.The Fresnel zone scope is vital to the excision reflection wave accurately, and by formula (19)-(21) as can be known, interface dip has material impact to the Fresnel zone scope.For addressing this problem, the scan method that the present invention uses the number percent increase and decrease obtains accurate reflector dip.Can the foundation of scanning be exactly reject reflection line-ups fully on the scattering wave imaging section.Whether this can by to the ratio deviation stacked section, observe continuous reflection line-ups and disappear to realize.

3. breakpoint scattering wave image enhancement

Add up by introducing for the weight coefficient of Fresnel zone and along angle, in theory, can not contained single big gun migration result of reflection wave; All single big gun migration result are cumulative, can not contained the scattering wave skew stacked section of reflection wave.Yet, when reception channel number weak when scattering wave or single shot record is less, after cumulative along angle remaining reflection wave and noise might with the same magnitude of scattering wave amplitude, so just be difficult to obtain clearly scattering wave imaging results.For this reason, need scattering wave is strengthened.By the skew road collection before cumulative along angle is done automatic gain by the road, realize balancing energy, can strengthen scattering wave at each seismic trace, because scattering wave has good same phasic property, the stack road concentrates scattering wave to be enhanced, and noise can't increase.Processing can strengthen the energy of scattering wave in the skew stacked section like this.The present invention introduces positive and negative to the angle of incident wave and reflection wave, by the expansion effective angle, promoted the actual effect that this enhancing is processed.

The present invention will utilize scattering wave that little breakpoint is carried out imaging, thereby better be familiar with minor fault and the breakpoint of underground structure.Be different from the medium small scale and change the scattering wave that causes, a key character of breakpoint scattering wave is that reversal of poles occurs by steady phase point the time scattering wave lineups.Utilize this characteristics, can by the imaging results of steady phase point both sides is subtracted each other, further give prominence to reflection wave; This subtraction can be suppressed the scattering wave imaging of other type and remaining reflection line-ups, plays an important role to promoting little breakpoint imaging.Its specific implementation process is as follows:

Concentrate the road of Fresnel zone both sides to superpose respectively in the two-dimensional migration road, obtain two imaging roads at each CDP point.And then can obtain two skew stacked sections that consisted of by Fresnel zone imaging road, left side and imaging road, Fresnel zone right side respectively.The amplitude of remembering respectively two skew stacked sections is A (x, T) and B (x, T), and wherein x is the horizontal coordinate that CDP is ordered, and T is time depth.Right | A (x, T) | and | B (x, T) | it is level and smooth to carry out large scale, obtains its low frequency component, is designated as E (x, T) and F (x, T).Only to the scattering wave skew stacked section of the independent imaging of little turn-off breakpoint then be

$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)$

ε is a stable factor a small amount of in the following formula.

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

Description of drawings

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

Fig. 2 is the isogram of the block migration velocity field, the South Sea that obtained by the actual seismic data.Numeral is the migration velocity value among the figure.Horizontal coordinate is pressed the whole coordinate definition of block among the figure, so do not start from scratch.

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

Fig. 4 is the angle dimension tangent plane of the CDP two-dimensional migration road collection of selecting 700 corresponding shot point coordinate 16.45km.As seen reflection line-ups is showed bending shape.Middle blank causes owing to lacking nearly offset distance seismologic record.

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

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

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

Fig. 7 is the skew stacked section that imaging road, Fresnel zone left side consists of.The conventional migration technique stacked section of comparison diagram 3, visible reflection line-ups is substantially disallowable.

Fig. 8 is the skew stacked section that imaging road, Fresnel zone right side consists of.The conventional migration technique stacked section of comparison diagram 3, visible reflection line-ups is substantially disallowable.

Fig. 9 is only to the scattering wave of the independent imaging of little turn-off breakpoint skew stacked section.Comparison diagram 7 and Fig. 8 as seen, remaining reflection line-ups and other interference are suppressed fully, the breakpoint of tomography is strengthened.

Figure 10 a and Figure 10 b have provided 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 that amplify Fig. 3 conventional migration technique result's part, and Figure 10 b is that amplify the part of the little breakpoint distribution plan of Fig. 9, and as seen contrast the position of little variation occurs at reflection line-ups from figure, and the breakpoint existence has just been indicated in the imaging of breakpoint scattering wave.

Figure 11 a and Figure 11 b have provided 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 that amplify Fig. 3 conventional migration technique result's part, and Figure 11 b is that amplify the part of the little breakpoint distribution plan of Fig. 9; The same with Figure 10, the position of little variation appears at reflection line-ups, and the existence of breakpoint has just been indicated in the imaging of breakpoint scattering wave.

Figure 12 has provided 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 that amplify Fig. 3 conventional migration technique result's part, and Figure 12 b is that amplify the part of the little breakpoint distribution plan of Fig. 9.The variation of reflection line-ups is not obvious among Figure 12 a, but the imaging of Figure 12 b breakpoint scattering wave explicitly points out the existence of breakpoint.This shows that this anti-brightly identify conventional reflection wave and be offset impalpable little turn-off breakpoint.

Embodiment

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

(1) records the seismic reflection signals of artificial epicenter excitation with the acquisition mode of single towing cable, be recorded on the tape.Concrete seismic reflection signals is to adopt the wall scroll towing cable collection, every towing cable 944 roads, track pitch 3.125m; At the front end of towing cable 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, Prestack seismic data is done conventional sound attenuation process; By making up a two dimensional cross-section along line direction and depth direction, seismic data is used the two-dimensional process method; For the part common midpoint, extract common midpoint gather, the common midpoint gather that extracts is made conventional NMO velocity picks up, the reciprocal value of resulting speed is done laterally average, with obtain laterally uniformly velocity field as initial offset speed.

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

(4) to CRP gather, utilize initial offset speed to do inverse dynamic correction, do again the speed after normal moveout correction obtains upgrading, the reciprocal value of the speed after each CDP point is upgraded is done space smoothing, and resulting level and smooth velocity field namely is the migration velocity field.Fig. 2 is the isogram of migration velocity field.Numeral is the migration velocity value among the figure.

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

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

(7) according to the migration velocity field, Prestack seismic data is done big gun territory pre-stack time migration, form the two-dimensional migration road collection that changes with shot point and angle at each CDP point.Fig. 4 is the skew road collection of the big gun record of the CDP point 700 shot point coordinate 16.45km of place, i.e. two-dimensional migration road collection angle dimension tangent plane.This big gun record has been expanded is the single shot record of 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 the shot point numbering, and the third dimension is angle, and angle changes to positive 60 degree from negative 60 degree; Prestack seismic data is pressed the common-shot-gather ordering; 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 the shot point coordinate common-shot-gather is divided into groups, with not on the same group common-shot-gather be placed on the different computing nodes of cluster computer; At each computing node, to whole seismic trace circulations; To each seismic trace, to the whole CDP dot cycles in the migration aperture; To each CDP point, determine initial imaging point by migration aperture, this imaging point below imaging point is calculated successively; First by shot point, geophone station and imaging point coordinate and migration velocity field, when determining to walk and the imaging weight coefficient, and then determine skew amplitude A by the seismologic record amplitude of this seismic trace; Then according to the offset distance of seismic trace, horizontal coordinate and the time depth of imaging point, pick up the start offset distance of the residual normal moveout of step 5 record and the excision that stretches; According to the time depth of offset distance and imaging point, judge whether skew amplitude A should be cut, if do not excise, the angle that then calculates shot point incident wave and reflecting interface normal is

${\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">cos</figure-callout>}}^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 the following formula, θ is the angle of shot point incident wave and reflecting interface normal, unit degree, and T is the time depth of imaging point, i.e. during outward journey, unit second, V _RmsBe imaging point (x ₀, the migration velocity of T) locating, unit meter per second, H are the offset distances of seismic trace, unit rice, p _sAnd p _rFor

$p_s=\frac{{(x_0-x_s)}^2}{V_{\mathrm{rms}}^2{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">T</figure-callout>}}^2},$ $p_r=\frac{{(x_0-x_g)}^2}{V_{\mathrm{rms}}^2{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">T</figure-callout>}}^2}$

X wherein ₀, x _gAnd x _sRespectively the horizontal coordinate of imaging point, geophone station and shot point, unit rice;

If Δ θ is the angular interval of predefined, unit degree is according to cos ²(k Δ θ-Δ θ 2)≤cos ²θ＜cos ²(k Δ θ-Δ θ/2) determines positive integer k; If the big gun under this seismic trace is numbered n, if the residual normal moveout that picks up is τ, its unit is second, works as x _g〉=x _s, then will be offset on (T+ τ, n, k Δ θ) element that amplitude A is added to the three-dimensional array of depositing migration result, namely on the positive-angle; Work as x _g＜x _s, then will be offset on (T+ τ, n ,-k Δ θ) element that amplitude A is added to the three-dimensional array of depositing migration result; Complete to the circulation of whole seismic traces, just obtain two-dimensional migration road collection with shot point and angle variation 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 of determining Fresnel zone, three-dimensional table is determined the angle domain Fresnel zone of each imaging point accordingly.

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

$a_1=\frac{\mathrm{\&rho;}({\mathrm{<figure-callout\; id="1"\; label="p"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">p</figure-callout>}}_1\sin \mathrm{\&phi;}\cos \mathrm{\&beta;}+{\mathrm{<figure-callout\; id="2"\; label="p"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">p</figure-callout>}}_2\sin \mathrm{\&beta;}\cos \mathrm{\&phi;})-\sin \mathrm{\&alpha;}(\cos \mathrm{\&beta;}+\cos \mathrm{\&phi;})}{\mathrm{\&rho;}\cos \mathrm{\&phi;}\cos \mathrm{\&beta;}({\mathrm{<figure-callout\; id="1"\; label="p"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">p</figure-callout>}}_1^2\cos \mathrm{\&phi;}-{\mathrm{<figure-callout\; id="2"\; label="p"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">p</figure-callout>}}_2^2\cos \mathrm{\&beta;})-2{\mathrm{<figure-callout\; id="1"\; label="p"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">p</figure-callout>}}_1\sin \mathrm{\&alpha;}\sin \mathrm{\&phi;}\cos \mathrm{\&beta;}+\left[2{\mathrm{<figure-callout\; id="3"\; label="p"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">p</figure-callout>}}_3\cos \mathrm{\&phi;}(\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}-\mathrm{\&rho;}{p}_2)\right]/{\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">cos</figure-callout>}}^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;}({\mathrm{<figure-callout\; id="1"\; label="p"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">p</figure-callout>}}_1^2\cos \mathrm{\&phi;}-{\mathrm{<figure-callout\; id="2"\; label="p"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">p</figure-callout>}}_2^2\cos \mathrm{\&beta;})-2{\mathrm{<figure-callout\; id="1"\; label="p"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">p</figure-callout>}}_1\sin \mathrm{\&alpha;}\sin \mathrm{\&phi;}\cos \mathrm{\&beta;}+\left[2{\mathrm{<figure-callout\; id="3"\; label="p"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">p</figure-callout>}}_3\cos \mathrm{\&phi;}(\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}-\mathrm{\&rho;}{p}_2)\right]/{\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">cos</figure-callout>}}^2\mathrm{\&beta;}}$

In upper two formulas

$\sin \mathrm{\&beta;}=\sin \left(2\mathrm{\&alpha;}\right)\sqrt{1/p^2-{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2\mathrm{\&phi;}}-\cos \left(2\mathrm{\&alpha;}\right)\sin \mathrm{\&phi;},$ $\cos \mathrm{\&beta;}=\sqrt{1-{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2\mathrm{\&beta;}},$ $\cos \mathrm{\&phi;}=\sqrt{1-{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2\mathrm{\&phi;}}$

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

${\mathrm{<figure-callout\; id="3"\; label="p"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">p</figure-callout>}}_3={\mathrm{<figure-callout\; id="1"\; label="p"\; filenames="HDA00002094878500041.png,HDA00002094878500042.png"\; state="\{\{state\}\}">p</figure-callout>}}_1\cos \mathrm{\&phi;}[\frac{\mathrm{\&rho;}\sin \left(2\mathrm{\&alpha;}\right)\sin \mathrm{\&phi;}\cos \mathrm{\&phi;}}{\sqrt{1-{\mathrm{\&rho;}}^2{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2}\mathrm{\&phi;}}+\cos \left(2\mathrm{\&alpha;}\right)\cos \mathrm{\&phi;}]$

The equal dimensionless of each variable in the 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 _sAnd imaging point (x ₀, T), make V _RmsBe the migration velocity at imaging point place, calculate

$\sin \mathrm{\&phi;}=\frac{x_0-x_s}{\sqrt{{(x_0-x_s)}^2+T^2{\mathrm{<figure-callout\; id="2"\; label="V\; rms"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">V</figure-callout>}}_{\mathrm{rms}}^{}}}$

X in the formula ₀And x _sUnit be rice, when T is outward journey, unit second, V _RmsUnit be meter per second.By sin φ, the ρ of imaging point place and sin α being rounded to determining deviation of three dimension directions showing, in table, obtain a according to the integer that obtains ₁And a ₂Make f _dBe the dominant frequency of seismic data, the 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;})}{{\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">cos</figure-callout>}}^2\mathrm{\&beta;}+{\mathrm{<figure-callout\; id="3"\; label="p"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">p</figure-callout>}}_3\sin \mathrm{\&beta;}{c}^{\&PlusMinus;}-0.5{\mathrm{<figure-callout\; id="3"\; label="p"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">p</figure-callout>}}_3^2{\left(c^{\&PlusMinus;}\right)}^2}$

C in the formula ^±Dimensionless,

Unit be rice; Be calculated as follows again shot point x _s, imaging point (x ₀, T-1/ (2f _d)) and acceptance point

Corresponding incident and reflection wave angle:

$\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">cos</figure-callout>}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 the formula

$p_s=\frac{{(x_0-x_s)}^2}{V_{\mathrm{rms}}^2{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">T</figure-callout>}}^2},$ $p_r^{\&PlusMinus;}=\frac{{(x_0-x_r^{\&PlusMinus;})}^2}{V_{\mathrm{rms}}^2{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="HDA00002094878500021.png,HDA00002094878500031.png"\; state="\{\{state\}\}">T</figure-callout>}}^2}$

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

The two-dimensional migration road collection of (9) each CDP being ordered is done balancing energy, determine first different time degree of depth upper angle territory Fresnel zone, then in angular deflection road collection excision Fresnel zone part, the remainder on Fresnel zone both sides is superposeed 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.

Two-dimensional migration road collection is done automatic gain by the road, realize balancing energy; If the horizontal coordinate that CDP is ordered is x ₀, shot point coordinate corresponding to road centralized shotpoint dimension is x _s, time depth is T, α is for the angle of skew lineups and surface level on this CDP point different time degree of depth, and take counter clockwise direction as positive-angle, V _Rms(T) be the migration velocity on 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 the used spacing of three-dimensional table in the step 8 sin φ, ρ and sin α are rounded, in table, pick up two real numbers according to the integer that obtains, and then by the angle domain Fresnel zone [φ on the method calculating different time degree of depth of step 8 ₁, φ ₂], its unit degree of being;

In the different time degree of depth, definition with the weight coefficient that angle changes 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

θ concentrates angle corresponding to angle dimension, unit degree for the skew road in the following formula.Angular-trace gather corresponding to each shot point concentrated in the two-dimensional migration road, be multiply by weight coefficient W (θ) at each time depth, just realized excision Fresnel zone part.With the two-dimensional migration road concentrate angle corresponding to whole shot points greater than with 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, by the scan method of number percent increase and decrease, whether eliminated reflection line-ups fully according to scattering wave skew stacked section, determine to be offset accurately the lineups inclination angle.

The angle of lineups and surface level carries out the number percent increase and decrease on the skew stacked section that step 6 is obtained, and to the universe angle value applying step 9 of each number percent increase and decrease, and with two imaging trace-stackings that each CDP point place obtains, forms stacked section; With the skew stacked section contrast that this section and step 5 obtain, if disappear at a certain imaging point reflection line-ups, the angle numerical value of this moment is exactly accurately inclination angles of these lineups.

(11) the accurate skew lineups inclination angle that utilizes step 10 to obtain, repeating step 9, obtain two skew stacked sections that consisted of by Fresnel zone imaging road, left side and imaging road, Fresnel zone right side respectively, this two section applied energies equilibrium is subtracted each other, only obtain the scattering wave skew stacked section to the independent imaging of little turn-off breakpoint.Fig. 7 and Fig. 8 are respectively the skew stacked sections that Fresnel zone imaging road, left side and imaging road, Fresnel zone right side consist of.Fig. 9 is only to the scattering wave of the independent imaging of little turn-off breakpoint skew stacked section.

Utilize accurately reflecting interface and horizontal plane angle, namely be offset the lineups inclination angle, repeating step 9, the imaging road on each CDP point Fresnel zone left side and right side is made up respectively, and two migration imaging sections that obtain remember that respectively its amplitude is A (x, T) and B (x, T), wherein x is the horizontal coordinate that CDP is ordered, and T is time depth; Right | A (x, T) | and | B (x, T) | it is level and smooth to carry out large scale, obtains its low frequency component, is designated as E (x, T) and F (x, T); Only to the scattering wave skew stacked section of the independent imaging of little turn-off breakpoint then be

$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)}$

ε is a stable factor a small amount of in the following formula.

(12) by software for display scattering wave is offset stacked section numerical value and is converted to the profile image that the little breakpoint of subsurface geological structure distributes, this profile image is that structure is along the tangent plane picture of line direction, this image can indicate the impalpable turn-off of conventional reflection wave migrated image less than minor fault and the breakpoint of 1/4 wavelength, be used for migration and the shutoff situation of understanding underground oil and gas, determine effective oil and gas reservoir.Figure 10, Figure 11 and Figure 12 have provided respectively the little breakpoint distributed image of Fig. 9 and the local contrast of Fig. 3 conventional migration technique stacked section.By Figure 10 and Figure 11 as seen, the position of little variation occurs at reflection line-ups, the existence of breakpoint has been indicated in the imaging of breakpoint scattering wave, and this shows that the present invention can identify breakpoint really.The variation of lineups is not obvious among Figure 12, and the existence of breakpoint has but been pointed out in the imaging of breakpoint scattering wave, and this shows this anti-bright impalpable little turn-off breakpoint of conventional reflection wave migrated image of identifying.The final sumbission of this processing is: the method can be indicated minor fault, the breakpoint of underground structure subtly, thereby can determine better effective oil and gas reservoir.

Claims (6)

1. the scattering wave prestack formation method of the little turn-off breakpoint of identification is characterized in that: may further comprise the steps: steps A, record the seismic reflection signals of artificial epicenter excitation with the acquisition mode of wall scroll towing cable or survey line, be recorded on the tape; Step B, read seismic signal from tape, form Prestack seismic data, Prestack seismic data is done conventional sound attenuation process; By making up a two dimensional cross-section along line direction and depth direction, seismic data is used the two-dimensional process method; For the part common midpoint, extract common midpoint gather, the common midpoint gather that extracts is made conventional NMO velocity picks up, the reciprocal value of resulting speed is done laterally average, with obtain laterally uniformly velocity field as initial offset speed; Step C, according to initial offset speed, carry out calculations of offset with conventional prestack time method, form CRP gather; Step D, to CRP gather, utilize initial offset speed to do inverse dynamic correction, do again the speed after normal moveout correction obtains upgrading, the reciprocal value of the speed after each CDP point is upgraded is done space smoothing, resulting level and smooth velocity field namely is the migration velocity field; Step e, according to the migration velocity field, Prestack seismic data is done conventional pre-stack time migration calculates, the CRP gather that obtains is done the moving school of high-order residue and the excision that stretches, the start offset distance of record residual normal moveout and the excision that stretches; Along offset distance stack CRP gather, obtain being offset stacked section; Step F, foundation skew stacked section and migration velocity field, the angle of lineups and surface level on definite skew stacked section is take counterclockwise as positive-angle; Step G, according to the migration velocity field, Prestack seismic data is done big gun territory pre-stack time migration, form the two-dimensional migration road collection that changes with shot point and angle at each CDP point; Step H, according to migration velocity field, shot point coordinate, imaging point coordinate, the imaging point place skew lineups angle with surface level, set up the three-dimensional table of definite Fresnel zone, three-dimensional table is determined the angle domain Fresnel zone of each imaging point accordingly; Step I, the two-dimensional migration road collection that each CDP is ordered are done balancing energy, determine first different time degree of depth upper angle territory Fresnel zone, then in angular deflection road collection excision Fresnel zone part, the remainder on Fresnel zone both sides is superposeed respectively, obtain two imaging roads at each CDP point; Step J, on the basis at the lineups inclination angle that step F obtains, by the scan method of number percent increase and decrease, according to whether having rejected reflection line-ups fully on the scattering wave skew stacked section, determine to be offset accurately the lineups inclination angle; Step K, the accurate skew lineups inclination angle that utilizes step J to obtain, repeating step I, obtain two skew stacked sections that consisted of by Fresnel zone imaging road, left side and imaging road, Fresnel zone right side respectively, this two section applied energies equilibrium is subtracted each other, only obtain the scattering wave skew stacked section to the independent imaging of little turn-off breakpoint; Step L, by software for display scattering wave is offset stacked section numerical value and is converted to the profile image that the little breakpoint of subsurface geological structure distributes, this profile image is that structure is along the tangent plane picture of line direction, this image can indicate the impalpable turn-off of conventional reflection wave migrated image less than minor fault and the breakpoint of 1/4 wavelength, be used for migration and the shutoff situation of understanding underground oil and gas, determine effective oil and gas reservoir.

2. a kind of scattering wave prestack formation method of identifying little turn-off breakpoint according to claim 1, it is characterized in that: in step G, described according to the migration velocity field, Prestack seismic data is done big gun territory pre-stack time migration, forming the two-dimensional migration road collection that changes with shot point and angle 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 the shot point numbering, the third dimension is angle, and angle changes to positive 60 degree from negative 60 degree; Prestack seismic data is pressed the common-shot-gather ordering; 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 the shot point coordinate common-shot-gather is divided into groups, with not on the same group common-shot-gather be placed on the different computing nodes of cluster computer; At each computing node, to whole seismic trace circulations; To each seismic trace, to the whole CDP dot cycles in the migration aperture; To each CDP point, determine initial imaging point by migration aperture, this imaging point below imaging point is calculated successively; First by shot point, geophone station and imaging point coordinate and migration velocity field, when determining to walk and the imaging weight coefficient, and then determine skew amplitude A by the seismologic record amplitude of this seismic trace; Then according to the offset distance of seismic trace, horizontal coordinate and the time depth of imaging point, pick up the start offset distance of the residual normal moveout of step e record and the excision that stretches; According to the time depth of offset distance and imaging point, judge whether skew amplitude A should be cut, if do not excise, the angle that then calculates shot point incident wave and reflecting interface normal is

In the following formula, θ is the angle of shot point incident wave and reflecting interface normal, unit degree, and T is the time depth of imaging point, i.e. during outward journey, unit second, V _RmsBe imaging point (x ₀, the migration velocity of T) locating, unit meter per second, H are the offset distances of seismic trace, unit rice, p _sAnd p _rFor

X wherein ₀, x _gAnd x _sRespectively the horizontal coordinate of imaging point, geophone station and shot point, unit rice;

If Δ θ is the angular interval of predefined, unit degree is according to cos ²(k Δ θ-Δ θ/2)≤cos ²θ＜cos ²(k Δ θ-Δ θ/2) determines positive integer k; If the big gun under this seismic trace is numbered n, if the residual normal moveout that picks up is τ, its unit is second, works as x _g〉=x _s, then will be offset on (T+ τ, n, k Δ θ) element that amplitude A is added to the three-dimensional array of depositing migration result, namely on the positive-angle; Work as x _g＜x _s, then will be offset on (T+ τ, n ,-k Δ θ) element that amplitude A is added to the three-dimensional array of depositing migration result; Complete to the circulation of whole seismic traces, just obtain two-dimensional migration road collection with shot point and angle variation 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 formation method of identifying little turn-off breakpoint according to claim 1, it is characterized in that: in step H, described according to the migration velocity field, the shot point coordinate, the imaging point coordinate, the angle of imaging point place skew lineups and surface level, set up the three-dimensional table of determining Fresnel zone, three-dimensional table is determined that angle domain Fresnel zone corresponding to each imaging point is achieved in that and is made that ρ is the ratio of imaging point place actual layer speed and migration velocity accordingly, α is the inclination angle of imaging point place reflecting interface, namely be offset the angle of lineups and surface level, take counterclockwise as positive-angle; Set up three-dimensional table take sin φ, ρ and sin α as three dimensions, deposit two dimensionless factors in the table, its value is calculated by following two formulas:

In upper two formulas

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

The equal dimensionless of each variable in the 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 _sAnd imaging point (x ₀, T), make V _RmsBe the migration velocity at imaging point place, calculate

X in the formula ₀And x _sUnit be rice, when T is outward journey, unit second, V _RmsUnit be meter per second.By sin φ, the ρ of imaging point place and sin α being rounded to determining deviation of three dimension directions showing, in table, obtain a according to the integer that obtains ₁And a ₂Make f _dBe the dominant frequency of seismic data, the unit hertz calculates

And then calculate two coordinate figures:

C ± dimensionless in the formula,

Unit be rice; Be calculated as follows again shot point x _s, imaging point (x ₀, T-1/ (2f _d)) and acceptance point

Corresponding incident and reflection wave angle:

In the formula

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

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

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

In the different time degree of depth, definition with the weight coefficient that angle changes is

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

φ ₁-5<θ<φ ₁

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

φ ₂＜θ<φ ₂+5

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

θ concentrates angle corresponding to angle dimension, unit degree for the skew road in the following formula.Angular-trace gather corresponding to each shot point concentrated in the two-dimensional migration road, be multiply by weight coefficient W (θ) at each time depth, just realized excision Fresnel zone part.With the two-dimensional migration road concentrate angle corresponding to whole shot points greater than with less than (φ ₁+ φ ₁The road of)/2 superposes respectively, just obtains two imaging roads at each CDP point.

5. a kind of scattering wave prestack formation method of 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, scan method by the number percent increase and decrease, according to whether having rejected reflection line-ups fully on the scattering wave skew stacked section, determine to be offset accurately the lineups inclination angle and be achieved in that on the skew stacked section that step F is obtained that the angle of lineups and surface level carries out number percent and increases and decreases, universe angle value applying step I to each number percent increase and decrease, and with two imaging trace-stackings that each CDP point place obtains, form stacked section; With the skew stacked section contrast that this section and step e obtain, if disappear at a certain imaging point reflection line-ups, the angle numerical value of this moment is exactly accurately inclination angles of these lineups.

6. a kind of scattering wave prestack formation method of identifying little turn-off breakpoint according to claim 1, it is characterized in that: in step K, the described accurate skew lineups inclination angle that utilizes step J to obtain, repeating step I, obtain two skew stacked sections that consisted of by Fresnel zone imaging road, left side and imaging road, Fresnel zone right side respectively, this two section applied energies equilibrium is subtracted each other, only obtaining scattering wave skew stacked section to the independent imaging of little turn-off breakpoint is achieved in that and utilizes accurately reflecting interface and horizontal plane angle, namely be offset the lineups inclination angle, repeating step I, the imaging road on each CDP point Fresnel zone left side and right side is made up respectively, and two migration imaging sections that obtain remember that respectively its amplitude is A (x, T) and B (x, T), wherein x is the horizontal coordinate that CDP is ordered, and T is time depth; Right | A (x, T) | and | B (x, T) | it is level and smooth to carry out large scale, obtains its low frequency component, is designated as E (x, T) and F (x, T); Only to the scattering wave skew stacked section of the independent imaging of little turn-off breakpoint then be

ε is a stable factor a small amount of in the following formula.

Images