CN101285894B - Heaved earth surface collected seismic data direct prestack time migration method - Google Patents

Abstract

The invention relates to a direct prestack time migration method for the seismic data acquired under the undulant ground surface. The method is applied to reflection seismic data processing during earthquake exploration and is a prestack migration imaging method aiming at the seismic data acquired under the undulant ground surface. The method can be directly applied to the seismic data which is acquired under the undulant ground surface and is acquired at shot points and demodulation points of different heights. Accordingly, the method can avoid errors caused by field static correction processing as well as correct the speed of a speed reducing belt in a migration process. The method can self-adaptively determine time-varying migration aperture and can correctly compensate geometric spreading effect of seismic waves during an imaging process so as to obtain an amplitude-preserved common reflection point gather. The core of the method is to use a one-way wave theory and a stable phase point principle to analyze so that the travel time and the amplitude of the seismic waves, which depend on the speed and thickness of the speed reducing belt, the heights of the shot points and the demodulation points, and stacking speed based on a floating datum plane, are obtained. The method has important application value for the explorations of oil, gas and mineral resources in areas with complex ground surfaces.

Description

The direct prestack time migration method of the seismic data that relief surface is gathered down

Technical field

The invention belongs to reflected seismic information processing technology field in the seismic prospecting, relate to the migration before stack imaging technique category in the seismic data processing procedure, is a kind of prestack time migration method at the seismic data of gathering under the relief surface.

Background technology

In the seismic data treatment scheme, the migration before stack imaging is the link of most critical, and pre-stack time migration is a kind of important method in the migration before stack imaging.Prestack time migration method can be to a class inclination angle, tomography complexity but the speed horizontal change is not the very violent better imaging of tectonic structure comparatively.Compare with prestack depth migration method, except that having higher counting yield, its main advantage is only to need to use stack velocity; Can obtain appropriate rate pattern by modes such as velocity sweepings simply like this, avoided a main difficulty using prestack depth migration method to face: speed is seen mould.Therefore, prestack time migration method has become the gordian technique of field of seismic exploration widespread use.

When seismic prospecting was carried out in surface conditions more complicated area, the collection of seismic data was to carry out on the face of land that rises and falls, and each shot point and geophone station be not on same elevation; And existing prestack time migration method is to ask for seismic travel-times with double square root equation, promptly uses pyramid face (pyramid) whilst on tour and superposes, and this requires shot point, geophone station is on same surface level.Solution to this problem has two kinds of disposal routes, and a kind of is to use static correcting method, promptly use static correction technology with all shot points and receiver static correction to same reference field, use the pre-stack time migration technology then; Second kind is also to be to use static correcting method, but will belong to the shot point of a CMP road collection and receiver static correction to the corresponding reference field of this CMP point, adopt " floating " rather than the same reference field of the whole district, Here it is floating datum prestack time migration method.The former is simple but static correction value is excessive, so error is also bigger; That use more at present is the latter.Although adopt the floating datum can be by reducing static correction value and mutual error counteracting reduces the static correction error, but still there is following several problem in this technology: the one, when high-velocity bed (old stratum) exposure, the static correction application conditions of vertical incidence and outgoing is false, continues to use static correcting method will bring bigger error; The 2nd, low velocity layer (LVL) is peeled off the error of replacing with speed and is difficult to offset fully; The 3rd, the static correction value known near-surface velocity model that places one's entire reliance upon can not be revised low velocity layer speed in imaging process.

The definite of migration aperture is important to pre-stack time migration, and the calculations of offset amount can be reduced in less aperture, but exists the risk that can not construct correct imaging to steep dip, and excessive aperture has been brought skew noise and bigger calculated amount again; Determine that better migration aperture remains the problem of a needs research.

Existing pre-stack time migration is based on double square root equation, and in fact its image-forming condition has used the dependent imaging condition of pre-stack depth migration, so imaging do not compensate the geometry diffusional effect of seismic wave propagation, does not promptly protect the width of cloth.Protect for obtaining common reflection point (CRP) the road collection of the width of cloth, utilize prestack information to carry out the technology of oil gas and fluid detection, should develop the formation method of protecting the width of cloth to serve prestack inversion etc.

Summary of the invention

The objective of the invention is: the prestack time migration method that the seismic data of gathering under a kind of relief surface is provided, can handle under the relief surface seismic data of gathering, result truly reflects the trend, breakpoint, turn-off size of underground structure and objectively to the reflection strength of seismic event.It does not promptly need to use the place static correcting method that shot point and geophone station are corrected on same elevation or the floating datum, can revise the speed of near surface low velocity layer again in migration process; This method has avoided in conjunction with the prestack time migration method of static correction technology because the error that static correction produces, and has also avoided obtaining in the static corrections processing difficulty of near-surface velocity field.

The technical solution used in the present invention is: the direct prestack time migration method of the seismic data that relief surface is gathered down, and concrete steps comprise:

(1) gathers the seismic data under the relief surface and read the seismic signal that each geophone station writes down with monolateral or bilateral acquisition mode, shot point is carried out migration processing with geophone station in the pre-stack seismic data of different elevations.

(2) by the information such as speed of known near surface, determine initial low velocity layer speed and low velocity layer time thickness, the low velocity layer bottom is defined as floating datum.

(3) adopt based on one way wave equation and the steady look-up method of putting principle mutually, try to achieve on the different elevations shot point and geophone station to the seimic travel time and the amplitude of imaging point, specifically: according to the low velocity layer speed v ₁, stack velocity v _RmsWith shot point or geophone station whilst on tour T to the horizontal range x of imaging point and imaging point to floating datum ₂With shot point or geophone station whilst on tour T to floating datum ₁Variation range, make up the table that seimic travel time and amplitude are tabled look-up and utilized in the algorithm, promptly by the T that equidistantly changes ₁v ₁, T ₂v _Rms, x and v ₁/ v _Rms, find the solution a=v _Rmsp _xWith amplitude A _sAnd they are filled out with T ₁v ₁, T ₂v _Rms, x and v ₁/ v _RmsBe in the four-dimension table of four parameters, its p _xBe the ray parameter that solves; By the low velocity layer thickness of each imaging point, speed, apart from the horizontal range of shot point and geophone station and based on the stack velocity of floating datum and the elevation of shot point and geophone station, can table look-up and try to achieve a and A _s, and then obtain shot point and geophone station seimic travel time and amplitude to imaging point.

The migration aperture that becomes when (4) determining according to the inclination maximum of intending the imaging structure, specifically: according to the information such as inclination maximum that are configured on the different depth (whilst on tour), become the initial whilst on tour table of migration aperture correspondence during structure,, can table look-up and try to achieve initial whilst on tour each imaging road.

(5) seismic trace is circulated, generate the imaging body of protecting the width of cloth by each seismic trace, specifically: to each imaging road, begin by initial whilst on tour, table look-up when trying to achieve shot point and geophone station and amplitude to the walking of each imaging point, use and protect the imaging amplitude that width of cloth image-forming condition obtains this point, common reflection point (CRP) road that then this amplitude is added to this imaging point correspondence is concentrated on the migration result of corresponding offset distance.

(6) in migration process, determine based on the stack velocity of floating datum and the speed of the initial low velocity layer of correction.

(7) with the skew of reforming of new stacking velocity field that obtains in the migration process and low velocity layer velocity field, all CRP road collection are evened up the CRP road collection that obtains protecting the width of cloth with the normal moveout correction method, collection stack in CRP road is obtained skew stacked section based on floating datum; Calculate the vertical whilst on tour T of floating datum and the level reference of selecting _f=T _c+ (h ₀-h _c)/v ₁, h wherein _cBe face of land elevation, T _cBe time thickness, the v of low velocity layer ₁Be revised low velocity layer speed, h ₀It is the elevation of the level reference selected; Use T _fRevise the skew stacked section, promptly press T _fEach seismic trace of mobile stacked section obtains final migration result.

(8) by conventional software for display migration result is converted to the profile image of underground reflective construct, profile image will react the trend, breakpoint, turn-off size of underground structure and to the reflection strength of seismic event.

Described employing based on one way wave equation and the look-up method of surely putting principle mutually try to achieve shot point and geophone station on the different elevations and be achieved in that to the seimic travel time of imaging point and amplitude and make v ₁Be the root-mean-square velocity of the low velocity layer more than the imaging point place floating datum, v _RmsBe the stack velocity of imaging point place based on floating datum, T ₂Be the vertical whilst on tour of imaging point to floating datum, x is shot point or the geophone station horizontal range to imaging point; If the face of land elevation at imaging point place is h _c, this place's low velocity layer time thickness be T _c, shot point or geophone station face of land elevation be h _s, T is arranged ₁=T _c+ (h _s-h _c)/v ₁With T ₁v ₁, T ₂v _Rms, x and v ₁/ v _RmsBe four parameters, set up four-dimensional table, deposit two numerical value in the table, a value is by equation:

$\frac{T_1{v}_1^2{p}_x}{\sqrt{1-v_1^2{p}_x^2}}+\frac{T_2{v}_{\mathrm{rms}}^2{p}_x}{\sqrt{1-v_{\mathrm{rms}}^2{p}_x^2}}+x=0$

The a=v that solves _Rmsp _x, another value is by a and b=v ₁/ v _RmsThe amplitude A of decision _s:

$A_s=\sqrt{\frac{2\mathrm{\&pi;}}{v_{\mathrm{rms}}c}},$ $c=\frac{T_1{b}^2}{(1-a^2{b}^2)\sqrt{1-a^2{\mathrm{<figure-callout\; id="2"\; label="b"\; filenames="S2008101141336E00042.png"\; state="\{\{state\}\}">b</figure-callout>}}^2}}+\frac{T_2}{(1-a^2)\sqrt{1-{\mathrm{<figure-callout\; id="2"\; label="a"\; filenames="S2008101141336E00042.png"\; state="\{\{state\}\}">a</figure-callout>}}^2}}$

To each imaging point, by T ₁, v ₁, T ₂, v _RmsCalculate four parameter T with x ₁v ₁, T ₂v _Rms, x and b=v ₁/ v _Rms, and by the table spacing round, can in table, check in numerical value a=v _Rmsp _xWith amplitude A _s, bring a into following formula and can get seimic travel time t _s:

$t_s={\mathrm{<figure-callout\; id="1"\; label="T"\; filenames="S2008101141336E00012.png"\; state="\{\{state\}\}">T</figure-callout>}}_1\sqrt{1-{\left(\mathrm{ab}\right)}^2}+{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="S2008101141336E00042.png"\; state="\{\{state\}\}">T</figure-callout>}}_2\sqrt{1-a^2}+\mathrm{ax}/v_{\mathrm{rms}}.$

The migration aperture that described inclination maximum according to plan imaging structure becomes when determining is achieved in that by each seismic trace and generates an elliptic imaging body, be configured in the initial whilst on tour that inclination maximum on the different whilst on tours determines the nonzero value in each imaging road according to intending imaging, the migration aperture that becomes during with the initial whilst on tour definition in imaging road; Initial whilst on tour T _bDetermine by following formula:

$T_b=\frac{2h}{v_{\mathrm{rms}}\left(\tan \right[2\mathrm{\&alpha;}\left(T_b\right)-\mathrm{\&theta;}]-\tan \mathrm{\&theta;})},$ $\tan \mathrm{\&theta;}=\frac{x-h}{T_b{v}_{\mathrm{rms}}}$

Wherein α is the inclination maximum on the different whilst on tours, v _RmsBe stack velocity, h is half offset distance, and x is the distance of imaging point apart from central point.Adopt the algorithm of tabling look-up in the practical application, depositing corresponding different offset distance, stack velocity and imaging point initial whilst on tour in the table, in table, pick up corresponding initial whilst on tour by the corresponding parameter in each imaging road apart from the imaging road of central point distance.

The imaging body that the width of cloth is protected in described generation is achieved in that order is t when having obtained shot point to the walking of imaging point _s, amplitude is A _s, geophone station is t during to the walking of imaging point _r, amplitude is A _r, make f (t) be the seismologic record that this big gun is cautious, imaging point place then, the imaging amplitude of protecting the width of cloth is:

$I\left(T_2\right)=\frac{A_r}{A_s}{f}_h(t_s+t_r)$

F in the formula _h(t) be half derivative of f (t) when two-dimensional case, promptly

With

Represent positive inverse-Fourier transform respectively; F during three-dimensional situation _h(t) be the first order derivative of f (t).

Describedly in migration process, determine to be achieved in that at first that based on the stack velocity of floating datum and the speed of revising the near surface low velocity layer seismic data is done preliminary velocity pick obtains initial, horizontal stacking velocity field uniformly, is offset based on this velocity field and initial low velocity layer speed; To being offset the CRP gather that generates, do inverse dynamic correction according to laterally uniform stacking velocity field, do the definite stack velocity of normal moveout correction again based on floating datum; The lineups of a selected high s/n ratio, increase and reduce initial low velocity layer speed based on new stacking velocity field with according to number percent, the neighborhood of these lineups is done local skew, choose the speed that makes lineups the most flat, this speed is done space smoothing handle, promptly obtain the low velocity layer velocity field of revising.

The direct prestack time migration method of the seismic data that relief surface of the present invention is gathered down can directly correctly be considered the influence to migration imaging of relief surface shape and near surface low velocity layer in migration process.

The direct prestack time migration method of the seismic data that relief surface of the present invention is gathered down, can in imaging process, revise low velocity layer speed, method is steady key to the selection of low velocity layer time thickness, does not require the high-velocity bed end face of its accurate corresponding actual configuration in bottom surface.

Relief surface of the present invention is the direct prestack time migration method of the seismic data of collection down, can correctly compensate the geometry diffusional effect of seismic event in imaging process, the CRP gather that obtains protecting the width of cloth.

The direct prestack time migration method of the seismic data that relief surface of the present invention is gathered down, the imaging of energy foundation plan is configured in the inclination maximum on the different depth, the migration aperture that becomes when determining adaptively.

The direct prestack time migration method of the seismic data that relief surface of the present invention is gathered down, but flexible Application is suitable for Parallel Implementation in all kinds of seismic channel sets.

Specific implementation principle of the present invention is as follows:

Core concept of the present invention is to analyze and the research pre-stack time migration from prestack depth migration method, regards pre-stack time migration as prestack depth migration method one approximate (the substituting the degree of depth with whilst on tour) of simplifying under the rate pattern.

At first based on the one way ripple theoretical with surely put principle mutually, obtain calculating and run point or geophone station to the seimic travel time of imaging point and the analytical expression of amplitude, this formula is based on (following discussion will with two-dimensional problems be example) of imaging point place based on the reduction of speed band root-mean-square velocity at the stack velocity of floating datum and floating datum top.

In wave number-frequency field, the degree of depth continuation of the seismic wave field of single shot point or geophone station can be expressed as:

$P(k_x,\mathrm{\&omega;},\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{t}_i)=F\left(\mathrm{\&omega;}\right)\exp (-j\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{t}_i\sqrt{{\mathrm{\&omega;}}^2-{\mathrm{\&alpha;}}_i^2{k}_x^2})---\left(1\right)$

Wherein F (ω) is the fourier-transform of the seismologic record of the focus of shot point or geophone station, t _iAnd α _iBe respectively the time thickness and the interlayer speed of each layer medium.Phase-shift phase in the formula (1) can be approximately

$\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{t}_i\sqrt{{\mathrm{\&omega;}}^2-{\mathrm{\&alpha;}}_i^2{k}_x^2}\&ap;{\mathrm{<figure-callout\; id="1"\; label="T"\; filenames="S2008101141336E00012.png"\; state="\{\{state\}\}">T</figure-callout>}}_1\sqrt{{\mathrm{\&omega;}}^2-v_1^2{k}_x^2}+{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="S2008101141336E00042.png"\; state="\{\{state\}\}">T</figure-callout>}}_2\sqrt{{\mathrm{\&omega;}}^2-v_{\mathrm{rms}}^2{k}_x^2}---\left(2\right)$

If establishing the m bed interface is floating datum, then ${\mathrm{<figure-callout\; id="1"\; label="T"\; filenames="S2008101141336E00012.png"\; state="\{\{state\}\}">T</figure-callout>}}_1=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{t}_i,$ ${\mathrm{<figure-callout\; id="1"\; label="v"\; filenames="S2008101141336E00012.png"\; state="\{\{state\}\}">v</figure-callout>}}_1^2=\frac{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i^2{t}_i}{{\mathrm{<figure-callout\; id="1"\; label="T"\; filenames="S2008101141336E00012.png"\; state="\{\{state\}\}">T</figure-callout>}}_1},$ ${\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="S2008101141336E00042.png"\; state="\{\{state\}\}">T</figure-callout>}}_2=\underset{i=m+1}{\overset{n}{\mathrm{\&Sigma;}}}{t}_i,$ $v_{\mathrm{<figure-callout\; id="2"\; label="rms"\; filenames="S2008101141336E00042.png"\; state="\{\{state\}\}">rms</figure-callout>}}^2=\frac{\underset{i=m+1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i^2{t}_i}{{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="S2008101141336E00042.png"\; state="\{\{state\}\}">T</figure-callout>}}_2}.$ Formula (2) shows that floating datum is the high-velocity bed end face of necessary corresponding actual configuration not; But if this reference field will effectively improve the approximation quality of formula (2) if accurately separate strong speed contrast interface (promptly separating low velocity layer and high-velocity bed end face).With the degree of depth continuation formula of approximate phase-shift phase substitution (1) formula of (2) formula and make space Fu Shi inverse transformation, the space-frequency domain wave field be

$P(x,\mathrm{\&omega;},{\mathrm{<figure-callout\; id="1"\; label="T"\; filenames="S2008101141336E00012.png"\; state="\{\{state\}\}">T</figure-callout>}}_1+T_2)=\frac{1}{2\mathrm{\&pi;}}\&Integral;F\left(\mathrm{\&omega;}\right)\exp [-\mathrm{j\&omega;}({\mathrm{<figure-callout\; id="1"\; label="T"\; filenames="S2008101141336E00012.png"\; state="\{\{state\}\}">T</figure-callout>}}_1\sqrt{1-{\mathrm{<figure-callout\; id="1"\; label="v"\; filenames="S2008101141336E00012.png"\; state="\{\{state\}\}">v</figure-callout>}}_1^2{(k_x/\mathrm{\&omega;})}^2}+{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="S2008101141336E00042.png"\; state="\{\{state\}\}">T</figure-callout>}}_2\sqrt{1-v_{\mathrm{rms}}^2{(k_x/\mathrm{\&omega;})}^2}-\frac{x{k}_x}{\mathrm{\&omega;}})]d{k}_x---\left(3\right)$

Make p _x=k _x/ ω, formula (3) can be exchanged into about p _xIntegration, promptly

$P(x,\mathrm{\&omega;},{\mathrm{<figure-callout\; id="1"\; label="T"\; filenames="S2008101141336E00012.png"\; state="\{\{state\}\}">T</figure-callout>}}_1+T_2)=\frac{\mathrm{\&omega;}}{2\mathrm{\&pi;}}\&Integral;F\left(\mathrm{\&omega;}\right)\exp [-\mathrm{j\&omega;}({\mathrm{<figure-callout\; id="1"\; label="T"\; filenames="S2008101141336E00012.png"\; state="\{\{state\}\}">T</figure-callout>}}_1\sqrt{1-{\mathrm{<figure-callout\; id="1"\; label="v"\; filenames="S2008101141336E00012.png"\; state="\{\{state\}\}">v</figure-callout>}}_1^2{p}_x^2}+{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="S2008101141336E00042.png"\; state="\{\{state\}\}">T</figure-callout>}}_2\sqrt{1-v_{\mathrm{rms}}^2{p}_x^2}-p_{x}x)]d{p}_x---\left(4\right)$

P in the formula _xTherefore be ray parameter,, can use (4) formula and surely put principle mutually and try to achieve asymptotic solution and be with frequency-independent

$P(x,\mathrm{\&omega;},{\mathrm{<figure-callout\; id="1"\; label="T"\; filenames="S2008101141336E00012.png"\; state="\{\{state\}\}">T</figure-callout>}}_1+T_2)=F\left(\mathrm{\&omega;}\right)\sqrt{\mathrm{\&omega;}}\exp (-j\frac{\mathrm{\&pi;}}{4})\sqrt{\frac{1}{2\mathrm{\&pi;}\left|{\mathrm{\&phi;}}^{\&prime;\&prime;}\left(p_x^0\right)\right|}}\exp (-\mathrm{j\&omega;\&phi;}\left(p_x^0\right))---\left(5\right)$

In the formula $\mathrm{\&phi;}\left(p_x\right)={\mathrm{<figure-callout\; id="1"\; label="T"\; filenames="S2008101141336E00012.png"\; state="\{\{state\}\}">T</figure-callout>}}_1\sqrt{1-{\mathrm{<figure-callout\; id="1"\; label="v"\; filenames="S2008101141336E00012.png"\; state="\{\{state\}\}">v</figure-callout>}}_1^2{p}_x^2}+{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="S2008101141336E00042.png"\; state="\{\{state\}\}">T</figure-callout>}}_2\sqrt{1-v_{\mathrm{rms}}^2{p}_x^2}-p_{x}x,{\mathrm{\&phi;}}^{\&prime;\&prime;}\left(p_x\right)$ Be its second derivative, p _x ⁰Be φ (p _x) zero point of first order derivative, φ (p _x ⁰) and $\sqrt{\frac{1}{2\mathrm{\&pi;}\left|{\mathrm{\&phi;}}^{\&prime;\&prime;}\left(p_x^0\right)\right|}}$ When promptly being the walking of seismic event and amplitude; And p _x ⁰Can try to achieve by following equation:

$\frac{T_1{v}_1^2{p}_x}{\sqrt{1-v_1^2{p}_x^2}}+\frac{T_2{v}_{\mathrm{rms}}^2{p}_x}{\sqrt{1-v_{\mathrm{rms}}^2{p}_x^2}}+x=0---\left(6\right)$

Formula (5) and (6) promptly are that shot point or geophone station are to the seimic travel time of imaging point and the analytical formula of amplitude, v in the formula ₁Be the root-mean-square velocity of this low velocity layer more than floating datum, v _RmsBe the stack velocity based on floating datum, T ₂Be the whilst on tour of imaging point apart from floating datum, T ₁Be face of land elevation h by the imaging point place _c, low velocity layer the time thickness T _c, shot point (or geophone station) face of land elevation h _s(or h _r) determine, T is arranged ₁=T _c+ (h _s-h _c)/v ₁

The present invention is based on direct rays hypothesis and ignore surface relief and determine migration aperture, below discussion will be example with two-dimensional problems.The present invention adopts the input channel imaging mode, and it is the elliptic imaging body of a focusing that each seismic trace will generate with shot point and geophone station, the time to become migration aperture be that initial whilst on tour with the nonzero value in each imaging road in the elliptic imaging body decides.If will be the structure imaging of α to the inclination angle, this imaging ellipse just should be more than or equal to α (semiellipse 90 degree with correspondence) at the maximum opening angle of this degree of depth; And will to incide this inclination angle be the structure of α and reflex to geophone station just from the seismic ray that shot point sends, and therefore determines the initial whilst on tour T of migration aperture correspondence according to following geometric relationship _h:

T _bv _rmstan[2α(T _b)-θ]-T _bv _rmstanθ＝2h (7)

T _bv _rmstanθ＝x-h

α is the inclination maximum that different depth (whilst on tour) is gone up the imaging structure in the formula, v _RmsBe the stack velocity of imaging point, θ is the incident angle of the seismic event of shot point outgoing, and 2 α-θ is the emergence angle of interface echo, and 2h is an offset distance, and x is the distance of imaging point apart from central point.Try to achieve T by (7) formula corresponding to different x _b, each imaging road is by T during imaging _bThe beginning imaging, the migration aperture that becomes when promptly realizing.

To each seismic trace since based on phase-shift method with surely put mutually principle tried to achieve shot point and geophone station to the amplitude of each imaging point when walking, utilize the deconvolution image-forming condition of Wave equation depth migration to obtain true amplitude imaging.If shot point and geophone station during to the walking of imaging point and amplitude be respectively t _s, A _sAnd t _r, A _r, be pulse between a period of time if establish focus, true amplitude imaging results (two-dimensional case) is arranged:

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

First three items in the integration can be regarded the fourier-transform of half derivative of this road earthquake record as, remembers that it is f in time domain _h(t), then (8) formula can be reduced to $I(x,T)=\frac{A_r}{A_s}{f}_h(t_s+t_r),$ Promptly in the seismologic record of doing after asking half derivative, pick up t _s+ t _rValue constantly is multiplied by the weight coefficient of being determined by amplitude then.Because real seismic record disperses, above-mentioned picking up can be realized by four point interpolations.

Beneficial effect of the present invention: this method can directly apply under the relief surface gather, shot point and the seismic data of geophone station at different elevations, generate pre-stack time migration road collection and image, this has not only save the processing links of place static correction, the error of also having been avoided static correction technology to bring.This method can be revised in migration process and reduce speed belt speed, and allows to reduce the high-velocity bed end face of the inaccurate corresponding actual configuration in speed belt bottom surface, and this has been avoided accurately obtaining in the disposal routes such as use static correction the difficulty of near-surface velocity model.This method can generate common reflection point (CRP) the road collection of protecting the width of cloth, serves oil gas and fluid detection technology such as prestack inversion better.The profile image that shows reflects the trend, breakpoint, turn-off size of underground structure and to the reflection strength of seismic event.This method has significant application value to oil gas, the mineral resources exploration of face of land complex area.

Description of drawings

Fig. 1 is the rifted-basin matter structural map with relief surface of gathering seismic data and carrying out pre-stack time migration, and the structure top layer does not have tangible weathering zone, has waste mantle about one deck maximum ga(u)ge 20m near surface, and speed is 1400m/s; The speed of other each several part is by (unit is m/s) shown in the numeral among the figure.We will carry out forward simulation by the full discrete values method of ACOUSTIC WAVE EQUATION based on this tectonic structure, generate the artillery simulators record, verify and explain the direct prestack time migration method of the seismic data of gathering under the relief surface with this.

Fig. 2 is the zero-offset section that extracts by in 151 big guns record of just drilling generation, and shallow-layer ground roll, the through reflection wave that involves the near surface weathering zone are disallowable.The reflection of horizontal interface in first reflection line-ups corresponding diagram 1 among the figure.White line is the floating datum of selecting.

Fig. 3 has provided at initial low velocity layer speed and the laterally uniformly initial stack velocity CDP=7500m CRP road collection of ordering after the match; As can be seen from Figure, there is the moving school of residue in lineups, also have the lineups of surface and interformational multiples correspondence.

Fig. 4 has provided the velocity sweeping result who after the match the CDP point of Fig. 3 is done the part skew in new stack velocity; Transverse axis is the speed chosen and the ratio of initial velocity among the figure; Can determine that from figure the best low velocity layer speed that this CDP is ordered should be 95% of initial low velocity layer speed.

Fig. 5 is the comparison of revised low velocity layer speed and initial even low velocity layer speed.

Fig. 6 has provided the CDP=7500m point identical with Fig. 3, the CRP road collection under stack velocity after the renewal and reduction of speed tape speed, and lineups are more straight as seen from the figure, focus on also better (the moving school of less residue amount can be revised in the stack flow process).

Fig. 7 is to be the skew stacked section in zero moment with the floating datum.

Fig. 8 is to be the skew stacked section of the correction of reference field with 264m place surface level.The skew lineups of Fig. 1 horizontal interface correspondence also become straight in Fig. 8, and this shows that revised low velocity layer speed is correct.Fig. 1 and Fig. 8 contrast as can be known, three larger tomographies, and two little structures are all by fine imaging, and section is also by imaging well.

Fig. 9 be Guizhou waterside town area actual acquisition seismic data migration result the part relatively.Last figure is the migration result that adopts conventional static correcting method, figure below is the migration result of direct prestack time migration method of the present invention, Fig. 9 cathetus section is the result of fault interpretation, and by fault interpretation as can be known, the inventive method has been portrayed trend and the breakpoint and the turn-off of underground structure well.

Embodiment

The direct prestack time migration method of the seismic data that relief surface is gathered down, rifted-basin tectonic structure at Fig. 1 with relief surface, the seismologic record that has synthesized the split shooting both sides collection of 151 big guns with numerical method, is example with geophone station in the pre-stack seismic data of different elevations with this group shot point, is specially following steps:

(1) the pre-stack seismic data that obtains at bilateral acquisition mode reads the seismologic record of each geophone station of each big gun record, and shot point is directly carried out migration processing with geophone station in the pre-stack seismic data of different elevations; Whole seismologic records are asked half derivative; Adopt parallel processing, by the offset distance scope seismologic record is assigned to each computer node, Fig. 2 is the seismologic record of zero-offset.

(2) determine initial low velocity layer speed and low velocity layer time thickness, the low velocity layer bottom is defined as floating datum.White line among Fig. 2 is the low velocity layer bottom that we define, and also be floating datum, and the time of white line on Fig. 2 is exactly low velocity layer time thickness.With horizontal even velocity (2313m/s) as initial low velocity layer speed.

(3) according to the low velocity layer speed v ₁, stack velocity v _RmsWith shot point or geophone station whilst on tour T to the horizontal range x of imaging point and imaging point to floating datum ₂With shot point or geophone station whilst on tour T to floating datum ₁Variation range, make up the table that seimic travel time and amplitude are tabled look-up and utilized in the algorithm, promptly by the T that equidistantly changes ₁v ₁, T ₂v _Rms, x and v ₁/ v _RmsThe combination of value, solving equation:

$\frac{T_1{\mathrm{<figure-callout\; id="1"\; label="v"\; filenames="S2008101141336E00012.png"\; state="\{\{state\}\}">v</figure-callout>}}_1^2{p}_x}{\sqrt{1-{\mathrm{<figure-callout\; id="1"\; label="v"\; filenames="S2008101141336E00012.png"\; state="\{\{state\}\}">v</figure-callout>}}_1^2{p}_x^2}}+\frac{T_2{v}_{\mathrm{rms}}^2{p}_x}{\sqrt{1-v_{\mathrm{rms}}^2{p}_x^2}}+x=0$

Obtain a=v _Rmsp _x, with a and b=v ₁/ v _RmsThe substitution following formula gets amplitude A _s:

$A_s=\sqrt{\frac{2\mathrm{\&pi;}}{v_{\mathrm{rms}}c}},$ $c=\frac{T_1{b}^2}{(1-a^2{b}^2)\sqrt{1-a^2{\mathrm{<figure-callout\; id="2"\; label="b"\; filenames="S2008101141336E00042.png"\; state="\{\{state\}\}">b</figure-callout>}}^2}}+\frac{T_2}{(1-a^2)\sqrt{1-a^2}}$

Then with a and A _sFill out with T ₁v ₁, T ₂v _Rms, x and v ₁/ v _RmsBe in the four-dimension table of four parameters.

(4) become the initial whilst on tour table of migration aperture correspondence when being configured in information architectures such as inclination maximum on the different depth (whilst on tour).Make α (T) be the inclination maximum on the different whilst on tours, v _RmsBe stack velocity, h is half offset distance, and x is the distance of imaging point apart from central point, by the v that equidistantly changes _Rms, h and x the combination of value, solving equation:

$T_b=\frac{2h}{v_{\mathrm{rms}}\left(\tan \right[2\mathrm{\&alpha;}\left(T_b\right)-\mathrm{\&theta;}]-\tan \mathrm{\&theta;})},$ $\tan \mathrm{\&theta;}=\frac{x-h}{T_b{v}_{\mathrm{rms}}}$

With the T that solves _hFill out with v _Rms, h and x be in the three-dimensional table of three parameters.In table, pick up corresponding initial whilst on tour T by the corresponding parameter in each imaging road during imaging _b

(5) to the seismic trace circulation, concrete steps comprise: 1) by the offset distance of seismic trace, the position of imaging point and stack velocity are tabled look-up and are determined the initial whilst on tour T in each road _b, if T _bGreater than the length of maximum imaging depth or seismologic record, show that imaging point has exceeded the imaging scope, need not carry out imaging.2) to each imaging point, by the low velocity layer speed v at imaging point top ₁With the time thickness T _c, the imaging point place face of land elevation h _c, shot point (or geophone station) face of land elevation h _s(or h _r) determine T ₁=T _c+ (h _s-h _c)/v ₁, make T ₂=T _b+ n Δ T (Δ T is the time depth sampling of the migrated section selected) is according to the reduction of speed tape speed v at imaging point place ₁With stack velocity v based on floating datum _Rms, and imaging point is apart from the horizontal range x of shot point and geophone station, algorithm obtains a=v by tabling look-up _Rmsp _xWith amplitude A _s, make b=v ₁/ v _Rms, bring a and b into following formula and can get seimic travel time t _s:

$t_s={\mathrm{<figure-callout\; id="1"\; label="T"\; filenames="S2008101141336E00012.png"\; state="\{\{state\}\}">T</figure-callout>}}_1\sqrt{1-{\left(\mathrm{ab}\right)}^2}+{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="S2008101141336E00042.png"\; state="\{\{state\}\}">T</figure-callout>}}_2\sqrt{1-a^2}+\mathrm{ax}/v_{\mathrm{rms}}$

In the time of can obtaining shot point and geophone station to the walking of imaging point like this and amplitude t _s, A _sAnd t _r, A _r(work as t _s+ t _rLength or T greater than seismologic record ₂Stop the increase of n during more than or equal to maximum imaging depth), use the image-forming condition of protecting the width of cloth and get imaging amplitude I (T ₂):

$I\left(T_2\right)=\frac{A_r}{A_s}{f}_h(t_s+t_r)$

F in the formula _h(t) be half derivative of this road earthquake record.3) according to the v that tables look-up and obtain _Rmsp _x, the seismic ray that calculating shot point and geophone station send and the intersection point of this imaging point place floating datum are tried to achieve the offset distance on this reference field in view of the above.4) according to the offset distance value, with imaging amplitude I (T ₂) the CRP road that is added to this imaging point correspondence concentrates on the result of corresponding offset distance.

(6) definite stack velocity and the initial low velocity layer speed of correction in migration process based on floating datum.Concrete steps comprise: 1) along the selected CMP road of coarse grid collection, by NMO technology and laterally average, determine a horizontal stacking velocity field uniformly, be offset based on this velocity field and initial low velocity layer speed that Fig. 3 is the CRP road collection that skew generates.2) whole CRP road collection is circulated, do inverse dynamic correction, again by the definite stacking velocity field of normal moveout correction based on floating datum according to horizontal uniformly initial stacking velocity field.3) along being selected to picture point than refined net, concentrate the lineups of a selected high s/n ratio by the CRP road of correspondence, increase and reduce initial low velocity layer speed based on new stacking velocity field with according to number percent, again the neighborhood of lineups is done local migration imaging, Fig. 4 is the result of the part skew of the collection part, a CRP road when increasing and reducing initial low velocity layer speed according to number percent.4) choosing that to make the most flat speed of lineups be revised low velocity layer speed, is 95% as determining number percent among Fig. 4, this speed is done space smoothing handle the low velocity layer velocity field that obtains revising.Consult Fig. 5.

(7) with the skew of reforming of new stacking velocity field and low velocity layer velocity field, Fig. 6 is the new CRP road collection that generates, and all CRP road collection are evened up common reflection point (CRP) the road collection that promptly obtains protecting the width of cloth with the normal moveout correction method; With the stack of CRP road collection, promptly obtain migration imaging section, as shown in Figure 7 based on floating datum; By face of land elevation h _c, low velocity layer the time thickness T _cWith revised low velocity layer speed v ₁, calculate floating datum and selected elevation at h ₀The vertical whilst on tour T of level reference _f=T _c+ (h ₀-h _c)/v ₁, each seismic trace of Fig. 7 is pressed T _fMove as the time, obtain migration result based on same level reference.

(8) migration result is converted to the profile image of underground reflective construct by software for display.Consult Fig. 8.

Claims (5)

1. the direct prestack time migration method of the seismic data gathered down of a relief surface, it is characterized in that adopting following steps: A) gather the seismic data under the relief surface and read the seismic signal that each geophone station writes down, shot point is carried out migration processing with geophone station in the pre-stack seismic data of different elevations with monolateral or bilateral acquisition mode; B) by known near surface information, determine initial low velocity layer speed and low velocity layer time thickness, the low velocity layer bottom is defined as floating datum; C) adopt based on one way wave equation and the steady look-up method of putting principle mutually and try to achieve on the different elevations shot point and geophone station to the seimic travel time and the amplitude of imaging point; The migration aperture that becomes when D) determining according to the inclination maximum of intending the imaging structure; E) to seismic trace circulation, generate the imaging body of protecting the width of cloth by each seismic trace, the amplitude of imaging body is added on the CRP gather by offset distance; F) in migration process, determine based on the stack velocity of floating datum and the speed of correction near surface low velocity layer; G) with the skew of reforming of final stacking velocity field that obtains in the migration process and low velocity layer velocity field, utilize the normal moveout correction method to eliminate the moving school of residue amount, migration result in the middle of obtaining, calculate the vertical whilst on tour of floating datum and the level reference of selecting, with migration result in the middle of this whilst on tour correction, obtain final migration result; H) by conventional software for display migration result is converted to the profile image of underground reflective construct, profile image reflects the trend, breakpoint, turn-off size of underground structure and to the reflection strength of seismic event.

2. the direct prestack time migration method of the seismic data that a kind of relief surface according to claim 1 is gathered down is characterized in that: described employing is tried to achieve shot point and geophone station on the different elevations based on one way wave equation and the steady look-up method of putting principle mutually and is achieved in that to the seimic travel time of imaging point and amplitude and makes v ₁Be the root-mean-square velocity of the low velocity layer more than the imaging point place floating datum, v _RmsBe the stack velocity of imaging point place based on floating datum, T ₂Be the vertical whilst on tour of imaging point to floating datum, x is shot point or the geophone station horizontal range to imaging point; T ₁Be shot point or geophone station whilst on tour to floating datum; If the face of land elevation at imaging point place is h _c, this place's low velocity layer time thickness be T _c, shot point or geophone station face of land elevation be h _s, T is arranged ₁=T _c+ (h _s-h _c)/v ₁With T ₁v ₁, T ₂v _Rms, x and v ₁/ v _RmsBe four parameters, set up four-dimensional table, deposit two numerical value in the table, a value is by equation:

$\frac{T_1{v}_1^2{p}_x}{\sqrt{1-v_1^2{p}_x^2}}+\frac{T_2{v}_{\mathrm{rms}}^2{p}_x}{\sqrt{1-v_{\mathrm{rms}}^2{p}_x^2}}+x=0$

The a=v that solves _Rmsp _x, p _xBe the ray parameter that solves, another value is by a and b=v ₁/ v _RmsThe amplitude A of decision _s:

$A_s=\sqrt{\frac{2\mathrm{\&pi;}}{v_{\mathrm{rms}}c}},$ $c=\frac{T_1{b}^2}{(1-a^2{b}^2)\sqrt{1-a^2{b}^2}}+\frac{T_2}{(1-a^2)\sqrt{1-a^2}}$

To each imaging point, by T ₁, v ₁, T ₂, v _RmsCalculate four parameter T with x ₁v ₁, T ₂v _Rms, x and b=v ₁/ v _Rms, and by the table spacing round, can in table, check in numerical value a=v _Rmsp _xWith amplitude A _s, bring a into following formula and can get seimic travel time t _s:

$t_s=T_1\sqrt{1-{\left(\mathrm{ab}\right)}^2}+T_2\sqrt{1-a^2}+\mathrm{ax}/v_{\mathrm{rms}}.$

3. the direct prestack time migration method of the seismic data that a kind of relief surface according to claim 1 is gathered down, it is characterized in that: the migration aperture that described inclination maximum according to plan imaging structure becomes when determining is achieved in that by each seismic trace and generates an elliptic imaging body, be configured in the initial whilst on tour that inclination maximum on the different whilst on tours determines the nonzero value in each imaging road according to intending imaging, the migration aperture that becomes during with the initial whilst on tour definition in imaging road; Initial whilst on tour T _bDetermine by following formula:

$T_b=\frac{2h}{v_{\mathrm{rms}}\left(\tan \right[2\mathrm{\&alpha;}\left(T_b\right)-\mathrm{\&theta;}]-\tan \mathrm{\&theta;})},$ $\tan \mathrm{\&theta;}=\frac{x-h}{T_b{v}_{\mathrm{rms}}}$

Wherein α is the inclination maximum on the different whilst on tours, v _RmsIt is stack velocity, h is half offset distance, x is the distance of imaging point apart from central point, adopt the algorithm of tabling look-up in the practical application, depositing corresponding different offset distance, stack velocity and imaging point initial whilst on tour in the table, in table, picking up corresponding initial whilst on tour by the corresponding parameter in each imaging road apart from the imaging road of central point distance.

4. the direct prestack time migration method of the seismic data that a kind of relief surface according to claim 1 is gathered down is characterized in that: the imaging body that the width of cloth is protected in described generation is achieved in that order is t when having obtained shot point to the walking of imaging point _s, amplitude is A _s, geophone station is t during to the walking of imaging point _r, amplitude is A _r, make f (t) be the seismologic record that this big gun is cautious, imaging point place then, the imaging amplitude of protecting the width of cloth is:

$I\left(T_2\right)=\frac{A_r}{A_s}{f}_h(t_s+t_r)$

F in the formula _h(t) be half derivative of f (t) when two-dimensional case, promptly

With

Represent positive inverse-Fourier transform respectively; F during three-dimensional situation _h(t) be the first order derivative of f (t).

5. the direct prestack time migration method of the seismic data that a kind of relief surface according to claim 1 is gathered down, it is characterized in that: describedly in migration process, determine to be achieved in that at first that based on the stack velocity of floating datum and the speed of revising the near surface low velocity layer seismic data is done preliminary velocity pick obtains initial, horizontal stacking velocity field uniformly, is offset based on this velocity field and initial low velocity layer speed; To being offset the CRP gather that generates, do inverse dynamic correction according to laterally uniform stacking velocity field, do the definite stack velocity of normal moveout correction again based on floating datum; The lineups of a selected high s/n ratio, increase and reduce initial low velocity layer speed based on new stacking velocity field with according to number percent, the neighborhood of these lineups is done local skew, choose the speed that makes lineups the most flat, this speed is done space smoothing handle, promptly obtain the low velocity layer velocity field of revising.

Images