CN101937100A - Pre-stack depth migration method - Google Patents

Abstract

The invention relates to the field of seismic prospecting, in particular to a pre-stack depth migration method. The method comprises the following steps of: performing migration velocity analysis on seismic data; calculating travel time, arc length, emergence angle and incidence angle by ray tracing; performing migration aperture calculation; and performing pre-stack depth migration by using a Kirchhoff integration vector migration formula. The method integrates migration velocity analysis, migration aperture selection, ray tracing and Kirchhoff integration formula, does not need wave field separation, and realizes multi-component simultaneous migration and accurate homing of converted wave. Based on the practical exploration data, the pre-stack depth migration method has good surface imaging effect and high imaging resolution, the continuity of the deep reflecting in-phase axis on a converted wave imaging section is improved, and the construction is clearer.

Description

A kind of prestack depth migration method

Technical field

The present invention relates to field of seismic exploration, concrete be meant a kind of prestack depth migration method.

Background technology

Skew is last step that geological data is handled, and also is most important, a crucial process, and is all the more so for the underground medium of complexity.

Because the non-horizontality of underground medium layer position, caused the earthquake stacked section can not reflect real layer configuration state, therefore need handle in a basic enterprising stepping line displacement of stack geological data, proofread and correct because there is the illusion time shift of the layer position reflection that is produced in stratigraphic dip, make the reflection playback of subterranean layer position, the true rolling shape and the occurrence on reflection stratum.

The migration processing method is divided into two classes, and a class is carried out on the basis of stack, is called post-stack migration; One class is handled in the stack line displacement of advancing, and is called migration before stack.The migration before stack precision will be higher than post-stack migration, and is better to the playback ability of complicated diffracted wave; But it is very big that the amount of calculation that needs is handled in migration before stack, obtained developing widely and using along with the raising of nearly computer technology during the last ten years.Wherein the difference according to offset method is handled in migration before stack, mainly is the difference of the time-space domain at method place, is divided into pre-stack time migration and pre-stack depth migration.Pre-stack depth migration is owing to the output section directly reflects the underground medium form of real space Depth Domain and obtains bigger concern in application.

Industry is considerably less about the pre-stack depth migration research of multi-component earthquake data at present, ripe relatively still time migration after stack of using.It is the tackling key problem focus and the difficult point of present field of seismic exploration that the migration before stack of multi-component earthquake data is handled, and current research uses with trial and also be confined to pre-stack time migration mostly, also is far from entering in the production practices application.And most of research also is confined on the basis of compressional wave skew, continue to use the thinking identical and carry out the migration processing of C ripple or equivalent transformed wave with pure compressional wave skew, space vector between multi-component seismic signal relation isolated fully come, carry out the migration processing of simple ripple, both lost contain in the wave field descend medium information galore, imaging effect is also relatively poor.

The Li Xiangyang of Univ Edinburgh UK professor leader's (the Hengchang Dai of research group, and Xiang-Yang Li.Effect of errors in the migration velocity model of PS-converted waves on traveltime accuracy in prestack Kirchhoff time migration in weak anisotropic media.Geophysics, 2008,73 (5): S195-S205) (Xiang-Yang Li, Jianxin Yuan.Converted-wave moveout and conversion-point equations in layered VTI media:Theory and applications.Journal of Applied Geophysics, 2003, the pre-stack time migration technology of many components transformed wave that 34:297-318) preliminary test is successful, preliminary test and application in some regional multicomponent seismic surveys, but the further tackling key problem of still needing is perfect.In addition, wear with Lee and the influence of migration result is carried out the discussion of numerical simulation and error at weak anisotropy medium converted shear wave pre-stack time migration medium velocity.They discussed respectively converted shear wave speed, vertical velocity than, effective velocity than the contribution of the precision of, anisotropic parameters to the whilst on tour error of calculation, improve for the imaging of time domain transformed wave and to have reference value preferably; And some thinkings can be used for reference by the pre-stack depth migration of transformed wave.

The wonderful moon of king's eighties in last century is at Kuo (Kuo; J.T.; and Dai; T.; 1984; Kirchhoff elastic wave migration for the case of non-coincident source and receiver:Geophysics; 49 (8); the thought of vector shift has been proposed on basis 1223-1238) the earliest; and elastic wave Kirchhoff integral formula (Qin Fuhao, Guo Yaxi, Wang Miaoyue, THE KIRCHHOFF ELASTIC WAVE MIGRATION METHOD, Chinese Journal of Geophysics have been released; 1988,31 (5): 577-587).This is many components while offset methods of the elastic wave under a kind of homogeneous isotropic medium hypothesis, during by the walking of PP, PS ripple during with two-component data inverse continuation to underground imaging point, the vibration mode on three components when obtaining each imaging point and just being excited.By the theogram of numerical simulation, they test this method, obtain result relatively preferably, have carried out preliminary exploration in this direction, but not further development so far.

This method has following two deficiencies: what 1) obtain is two total component vibrational energies, does not realize wave field separation; 2) migration result that obtains of this method is the vibrations that the energy transmission of a certain shot point excites to for subsurface reflective boundary the time.Vibration vector when different shot points excite for underground same reflection spot is inequality, therefore can not be used for the stack of many big gun records.And their research is a kind of approximate pre-stack depth migration, promptly in migration process, still carry out in time domain, when just carrying out afterwards-dark conversion converting the migrated section of output to depth section.And the calculating when walking they only rest on single interface direct rays and follow the trail of this aspect, do not consider the difference of migration aperture for different Depth Imaging yet, the acquisition of migration velocity is simple relatively, just simply use the conversion that stack velocity is carried out interval velocity then, thereby effect is unsatisfactory.

Publication number is a kind of technology in the invention patent protection migration velocity analysis field of US 2009/0257308A1, and this technology is come continuous renewal speed model based on the normal moveout correction error, thereby makes final rate pattern approach the true velocity of underground medium.This patent is based on the velocity analysis technology of residue curvature analysis method, and this technology promptly requires the residue whilst on tour to can be approximated to be hyperbolic curve or para-curve based on little geophone offset hypothesis and flat seam position medium hypothesis.And under the violent or big geophone offset situation of speed horizontal change, this technology institute based on algorithm be similar to inappropriate; And this patented technology is only carried out Treatment Analysis at conventional compressional wave data and velocity of longitudinal wave.

(QiXue Dong, Jia Guanghua, Tong Zhaoqi such as QiXue Dong, Kirchhoff pre-stack time migration processing technique and application, physical prospecting geochemical exploration computing technique, 2009,31 (2): 126-130) research is pre-stack time migration technology and application thereof at conventional compressional wave Data Processing, and this technology has been ripe use, business-like production technology in the oil-gas seismic exploration field.

(Zhang Meng, Meng Xiangbin, Kuang Bin etc. such as Zhang Meng, integral method pre-stack depth migration technology is in the application in BS6 area, the exploration geophysics progress, 2009,32 (1): 48-55) introduce to be the pre-stack depth migration technology of conventional longitudinal wave earthquake data of ripe, commercial applications and software thereof carry out Processing Test and application to the geological data of Shengli Oil Field block, and carry out the contrast of imaging effect, thereby the advantage and the effect of explanation pre-stack depth migration with the post-stack migration (time) of conventional compressional wave.

(Ye Yueming, Li Zhenchun, Tong Zhaoqi such as Ye Yueming, based on guarantor's width of cloth pre-stack depth migration of stablizing image-forming condition, geophysical prospecting for oil, 2009,44 (1): 28-32) be primarily aimed at the problem that damp type deconvolution processing exists in the conventional compressional wave pre-stack depth migration processing, caused calculating wild effect improved the computational stability of migration processing by Gauss's smoothing processing when denominator was gone to zero, and reached the purpose that the part pressure is made an uproar simultaneously.

Above-mentioned relevant technologies at geological data all be simple component, conventional compressional wave data, none relates to the processing of multi-component earthquake data, related technology and theoretical method are not suitable for the processing of three-component seismic data yet.

Publication number is that the patent of invention of CN 101598805A is at the compressional wave and the difference in converted shear wave travel-time, by the method for speed than scanning and tracing of horizons, to demarcate, pick up with the different reflected phase will of converted shear wave section at compressional wave with one deck position, thereby the p-and s-wave velocity ratio of each corresponding phase is decided, and then realize the converted shear wave section according to the longitudinal wave propagation Time Compression, make on identical longitudinal wave propagation time scale, make things convenient for the layer position contrast of compressional wave and converted shear wave.The content of its protection does not relate to the imaging technique of three-component seismic data.

Publication number is because the problem that mass data processing is consuming time, counting yield is low during the patent of invention of CN 1797031A is handled at the pre-stack depth migration of conventional compressional wave data, coupling and combination by the same phase big gun, realize that many big guns merge a big gun, thereby reduce the big gun collection number of migration processing, on the basis of not damaging offset effect, realize the raising of migration processing efficient.The content of its protection does not relate to the imaging technique of three-component seismic data.

Publication number is that the content of the invention patent protection of CN 101419292A is the geological data at three-component acquisition, and the treatment scheme by prestack pre-service, static correction, velocity analysis, normal moveout correction, the such routine of stack is processed into picture to converted shear wave.The content of this patent protection generally is called conventional poststack disposal route again in the three-component seismic data process field, only relates to the time migration technology.

Publication number is that the content of the invention patent protection of US 2004/0117123A1 is the pre-stack time migration processing technique at the pure compressional wave data of routine, and this technology realizes at frequency-wavenumber domain.This patent is that geological data is transformed to frequency-wavenumber domain, and the linear fit that figure carries out phase shift to the ray parameter of given space any point when walking according to travelling then makes the violent complex dielectrics imaging problem of horizontal change of pre-stack time migration algorithm speed-adaptive.

Publication number is that the content of the invention patent protection of US 2010/0054082A1 does not relate to concrete migration algorithm, and is more irrelevant with many components pre-stack depth migration.All relate to the simulation of wave field in any migration process, it is not that calculating is finished when walking by ray tracing, but realizes by the finite-difference modeling of ACOUSTIC WAVE EQUATION.This patent is primarily aimed at that the prestack reverse-time migration needs the wave field of storage space any point in handling and a large amount of consumption problem of causing internal memory and hard drive space in the migration processing, and make and be offset inefficiency, a kind of technological improvement measure of saving calculator memory and hard drive space is proposed for this reason, promptly in the simulation of just drilling wave field, only preserve the wave field value on the impedance interface point, other spatial point do not participate in calculating, thereby improved the speed of migration processing greatly, reduced the consumption of computational resource.

Pre-stack depth migration is the effective tool to the violent complex structure imaging of big stratigraphic dip and speed horizontal change.But pre-stack depth migration relies on more accurate rate pattern.Conventional speed modeling method can't reach enough precision.Simultaneously, because pre-stack depth migration is very responsive to rate pattern, so it is as a kind of strong velocity analysis instrument, develop and two kinds of migration velocity analysis methods: deep focus analysis (DFA, depth-focusing analysis) and the residue curvature analysis (RCA, residual-curvature analysis).Whether accurately carrying out migration velocity analysis just must solve two problems: 1) set up evaluation speed standard; 2) how to carry out Velocity Updating.

The principle of deep focus analysis is according to stack power measuring speed error.When the skew degree of depth was identical with the depth of focus, speed reached requirement.The residue curvature analysis is according to residue travel-time difference measuring speed error.When the imaging depth of different geophone offsets was identical, speed reached requirement.

Because velocity analysis is normally based on a coarse primary iteration speed, therefore the velocity analysis method of a stalwartness also will guarantee the convergence of speed when reducing iterations.Traditional velocity analysis method is normally based on following three hypothesis: 1. stratiform uniform dielectric; 2. little geophone offset; 3. acline.

(MacKay S.and Abma R.Imaging and velocity estimation with depth-focusing analysis.Geophysics.1992, DFA method 57:1608-1622) has been used whole above three and has been similar to for MacKay and Abma.

In RCA, the offset method difference that is adopted during velocity analysis, based on hypothesis also different: 1. shot record migration or common receiver skew: (Al-Yahya is analysis by iterative profile migration K.1989.Velocity for Al-Yahay, Geophysics, 1989,54:718-729) method of Ti Chuing is based on above three hypothesis; Lee and Zhang (Lee W.and Zhang L.Residual shot profile migration, Geophysics, 1992, method 57:815-822) is extended to small inclination with the 3rd hypothesis.2. geophone offset skew altogether: (Deregowsiki S.M.Common-offset migration and velocity analysis.First Break, 1990,8 (6): method 225-234) has been used preceding two hypothesis to Deregowski.

Above method is all based on little geophone offset hypothesis, and under little geophone offset situation, the residue whilst on tour can be approximated to be hyperbolic curve or para-curve.When the speed horizontal change was violent, it was inappropriate using this approximate.Therefore, Liu (Liu Z.and Bleistein N.Migration velocity analysisi:Theory and an iterative algorithm, Geophysics, 1995,60:142-153) at Lafond and Levander (Lafond C.F.and Levander A.R.Migration moveout analysis and depth focusing, Geophysics, 1993, a kind of new RCA method has been proposed on basis 58:91-100)---based on the common imaging point velocity analysis method of deep focus, this method is applicable to any geophone offset, inclination angle and velocity variations.

Summary of the invention

Technical matters to be solved by this invention is at multi-component earthquake data, and a kind of three-component prestack depth migration method based on the vector wave field is provided, and does not need to carry out wave field separation, and realizes many component skews, realizes the accurate playback of transformed wave.

The technical scheme that the present invention solves the problems of the technologies described above is as follows: a kind of prestack depth migration method may further comprise the steps:

Step 1, geological data is carried out migration velocity analysis;

Step 2, when calculating away by ray tracing, arc length, emergence angle, incident angle;

Step 3, carry out migration aperture and calculate;

Step 4, utilize Kirchhoff integral method vector shift formula to carry out pre-stack depth migration.

The invention has the beneficial effects as follows: this method with migration velocity analysis, migration aperture choose, ray tracing, Kirchhoff integral formula be incorporated into one, do not need to carry out wave field separation, and realize many components skews, realizes the accurate playback of transformed wave.

On the basis of technique scheme, the present invention can also do following improvement.

Further, migration velocity analysis described in the step 1 adopts the common imaging point velocity analysis method based on deep focus.

Adopting the beneficial effect of above-mentioned further scheme to be, is to revise to form on residue curvature analysis method based on the common imaging point velocity analysis method of deep focus, and its fast convergence rate is applicable to various rate patterns and geophone offset.

Further, described step 1 comprises the steps:

Step a1, initial geologic model, and enter step a2;

Step a2, set up initial velocity model, and enter step a3;

Step a3, carry out the PSDM (pre-stack depth migration, pre-stack depth migration) of based target layer position bottom boundary, and enter step a4;

Step a4, carry out migration velocity analysis, and enter step a5;

Step a5, judged whether to set up complete rate pattern, if for "Yes" then enter step 2, if for "No" then enter step a6;

Step a6, carry out the PSDM of based target layer position bottom boundary, and enter step a7;

Step a7, pick up the destination layer bottom boundary, and enter step a8;

Step a8, renewal speed model move to down one deck position with destination layer, and enter step a3.

Wherein, complete rate pattern is described in the step a5: from shallow complete in speed, density, the layer position buried depth seismogeology model as parameter to dark each layer position.

Further, in the step 2, ray casting when adopting direct rays back tracking method or employing to walk based on the minimum of Fermat principle.

Adopt the beneficial effect of above-mentioned further scheme to be, adopt direct rays back tracking method computing velocity fast, ray casting can be avoided the ray blind area when employing was walked based on the minimum of Fermat principle.

Further, employing is applicable to that the elastic wave Kirchhoff integral method vector shift formula of heterogeneous anisotropic media carries out pre-stack depth migration in the step 4.

Adopt the beneficial effect of above-mentioned further scheme to be, be applicable to the elastic wave Kirchhoff integral method vector shift formula of heterogeneous anisotropic media, have the characteristic of wave field separation, separate PP ripple and PS ripple when can accomplish to be offset; Migration result can be exported PP wave profile and PS wave profile respectively, is beneficial to AVO (Amplitude Versus Offset, amplitude is with the variation of the offset distance) feature of analyzing PP ripple and PS ripple respectively; Use the Green function of elastic oscillation equation, improved guarantor's width of cloth of algorithm; Be applicable to heterogeneous anisotropic media.

The present invention is at people such as Wang Miaoyue (Qin Fuhao, Guo Yaxi, Wang Miaoyue; THE KIRCHHOFF ELASTIC WAVE MIGRATION METHOD, Chinese Journal of Geophysics, 1988; 31 (5): 577-587) release new offset equation on Yan Jiu the basis, and improve at its existing two deficiencies.In migration process, directly export PP ripple and PS ripple, when carrying out many big guns stacks, realized the effect of relevant reinforcement, obtain the section of PP and PS ripple.

Description of drawings

Fig. 1 is the schematic diagram one of prestack depth migration method of the present invention;

Fig. 2 is the schematic diagram two of prestack depth migration method of the present invention;

Fig. 3 is the schematic diagram three of prestack depth migration method of the present invention;

Fig. 4 is the process flow diagram of prestack depth migration method of the present invention;

Fig. 5 is the realization flow of migration velocity analysis among the present invention;

Fig. 6 A is for adopting based on the PP ripple design sketch after the time migration after stack method of CCP (Common Converted Point, transfer point altogether) road collection is offset 2D3C geological data of grand celebration altogether;

Fig. 6 B is for adopting based on the PS ripple design sketch after the time migration after stack method of CCP road collection is offset 2D3C geological data of grand celebration altogether;

Fig. 7 A is the PP ripple design sketch after the migration processing that adopts the inventive method a 2D3C geological data of grand celebration is carried out;

Fig. 7 B is the PS ripple design sketch after the migration processing that adopts the inventive method a 2D3C geological data of grand celebration is carried out;

Fig. 8 is the simulated data model of prestack depth migration method of the present invention;

Fig. 9 A is the geological data Z component record diagram of the data model of Fig. 8;

Fig. 9 B is the geological data X component record diagram of the data model of Fig. 8;

The compressional wave migrated depth section result of Figure 10 for adopting prestack depth migration method of the present invention to obtain in the data model of Fig. 8;

The converted shear wave migrated depth section result of Figure 11 for adopting prestack depth migration method of the present invention to obtain in the data model of Fig. 8;

The compressional wave migrated depth section result of Figure 12 for adopting the ACOUSTIC WAVE EQUATION skew in the data model of Fig. 8 and directly obtaining with the inclination angle perforate;

The converted shear wave migrated depth section result of Figure 13 for adopting the ACOUSTIC WAVE EQUATION skew in the data model of Fig. 8 and directly obtaining with the inclination angle perforate;

The compressional wave migrated depth section result of Figure 14 for adopting ACOUSTIC WAVE EQUATION skew and constant aperture to obtain in the data model of Fig. 8;

The converted shear wave migrated depth section result of Figure 15 for adopting ACOUSTIC WAVE EQUATION skew and constant aperture to obtain in the data model of Fig. 8.

Embodiment

Below in conjunction with accompanying drawing principle of the present invention and feature are described, institute gives an actual example and only is used to explain the present invention, is not to be used to limit scope of the present invention.

The method principle

As shown in Figure 1, the trivector wave field extrapolation equation of homogeneous isotropic medium is:

$u_n(\stackrel{\&RightArrow;}{x},t)=\frac{1}{4\mathrm{\&pi;}}\underset{-\&infin;}{\overset{\&infin;}{\text{\&Integral;}}}\underset{-\&infin;}{\overset{\&infin;}{\&Integral;}}\left\{\begin{array}{c}[(1-{2v}^2){\mathrm{\&delta;}}_{m3}{\hat{r}}_n+{2v}^2{\hat{r}}_m{\hat{r}}_n{\hat{r}}_3]\&times;{\stackrel{\&CenterDot;}{u}}_m({\stackrel{\&RightArrow;}{x}}^{\&prime;},t+t_P)/\left({\mathrm{rV}}_P\right)\\ +({\mathrm{\&delta;}}_{\mathrm{mn}}{\hat{r}}_3+{\mathrm{\&delta;}}_{n3}{\hat{r}}_m-2{\hat{r}}_m{\hat{r}}_n{\hat{r}}_3)\&times;{\stackrel{\&CenterDot;}{u}}_m({\stackrel{\&RightArrow;}{x}}^{\&prime;},t+t_S)/\left({\mathrm{rV}}_S\right)\\ +\left({\hat{r}}_m{\hat{r}}_n\right)T_{3m}({\stackrel{\&RightArrow;}{x}}^{\&prime;},t+t_P)/\left({\mathrm{\&rho;rV}}_P^2\right)\\ +({\mathrm{\&delta;}}_{\mathrm{mn}}-{\hat{r}}_m{\hat{r}}_n)T_{3m}({\stackrel{\&RightArrow;}{x}}^{\&prime;},t+t_S)/\left({\mathrm{\&rho;rV}}_S^2\right)\end{array}\right\}{\mathrm{<figure-callout\; id="1"\; label="dx"\; filenames="HSA00000232213300041.png,HSA00000232213300042.png"\; state="\{\{state\}\}">dx</figure-callout>}}_1^{\&prime;}{\mathrm{dx}}_2^{\&prime;}---\left(1\right)$

Wherein

The expression particle displacement, unit is m/s;

The expression particle velocity, the m of unit; Expression stress, unit is N/m ²M, n=1,2,3 represent x respectively, y, z coordinate, the m of unit;

Be the unit radius vector,

Represent the projection of vector of unit length on m, n coordinate of ray emergence angle respectively; T=t (x|x ') expression focus is during to the walking of underground diffraction point, and unit is s; t _p, t _sRepresent respectively diffraction put the PP ripple of acceptance point, when the PS ripple is walked, unit is s; R represents that diffraction puts the distance of acceptance point, and unit is m; V _p, V _sRepresent longitudinal and transverse wave velocity respectively, unit is m/s; V=V _s/ V _pExpression p-and s-wave velocity ratio.

Because Free Surface boundary condition, face of land stress are zero, can cast out, have for back two of equation (1)

$u_n(\stackrel{\&RightArrow;}{x},t)=\frac{1}{4\mathrm{\&pi;}}\underset{-\&infin;}{\overset{\&infin;}{\&Integral;}}\underset{-\&infin;}{\overset{\&infin;}{\&Integral;}}\left\{\begin{array}{c}[(1-{2v}^2){\mathrm{\&delta;}}_{m3}{\hat{r}}_n+{2v}^2{\hat{r}}_m{\hat{r}}_n{\hat{r}}_3]\&times;{\stackrel{\&CenterDot;}{u}}_m({\stackrel{\&RightArrow;}{x}}^{\&prime;},t+t_P)/\left({\mathrm{rV}}_P\right)\\ +({\mathrm{\&delta;}}_{\mathrm{mn}}{\hat{r}}_3+{\mathrm{\&delta;}}_{n3}{\hat{r}}_m-2{\hat{r}}_m{\hat{r}}_n{\hat{r}}_3)\&times;{\stackrel{\&CenterDot;}{u}}_m({\stackrel{\&RightArrow;}{x}}^{\&prime;},t+t_S)/\left({\mathrm{rV}}_S\right)\end{array}\right\}{\mathrm{<figure-callout\; id="1"\; label="dx"\; filenames="HSA00000232213300041.png,HSA00000232213300042.png"\; state="\{\{state\}\}">dx</figure-callout>}}_1^{\&prime;}{\mathrm{dx}}_2^{\&prime;}---\left(2\right)$

This equation is launched, the energy on three components is projected to respectively on PP component (3) formula and PS component (4) formula, as shown in Figure 2

u _PP＝u ₁cosθsinφ+u ₂sinθsinφ+u ₃cosφ(3)

u _PS＝u ₁cosθcosφ+u ₂sinθcosφ+u ₃sinφ(4)

To bring (3) formula, (4) formula into after the expansion of (2) formula,

$u_{\mathrm{PP}}(x,y,z,t)=\frac{1}{4\mathrm{\&pi;}}\underset{-\&infin;}{\overset{\&infin;}{\&Integral;}}\underset{-\&infin;}{\overset{\&infin;}{\&Integral;}}\left[\begin{array}{c}\frac{{2v}^2{\hat{r}}_x{\hat{r}}_z}{{\mathrm{rV}}_P}\&CenterDot;{\stackrel{\&CenterDot;}{u}}_x(x^{\&prime;},y^{\&prime;},0,t+t_P)+\frac{{2v}^2{\hat{r}}_y{\hat{r}}_z}{{\mathrm{rV}}_P}\&CenterDot;{\stackrel{\&CenterDot;}{u}}_y(x^{\&prime;},y^{\&prime;},0,t+t_P)\\ +\frac{1-{2v}^2+{2v}^2{\hat{r}}_z{\hat{r}}_z}{{\mathrm{rV}}_P}\&CenterDot;{\stackrel{\&CenterDot;}{u}}_z(x^{\&prime;},y^{\&prime;},0,t+t_P)\end{array}\right]{\mathrm{dx}}^{\&prime;}{\mathrm{dy}}^{\&prime;}---\left(5\right)$

$u_{\mathrm{PS}}(x,y,z,t)=\frac{1}{4\mathrm{\&pi;}}\underset{-\&infin;}{\overset{\&infin;}{\&Integral;}}\underset{-\&infin;}{\overset{\&infin;}{\&Integral;}}\left[\begin{array}{c}\frac{1-2{\hat{r}}_z{\hat{r}}_z}{{\mathrm{rV}}_S}\&CenterDot;{\stackrel{\&CenterDot;}{u}}_x(x^{\&prime;},y^{\&prime;},0,t+t_S)+\frac{1-2{\hat{r}}_z{\hat{r}}_z}{{\mathrm{rV}}_S}\&CenterDot;{\stackrel{\&CenterDot;}{u}}_y(x^{\&prime;},y^{\&prime;},0,t+t_S)\\ +\frac{2({\hat{r}}_x{\hat{r}}_z+{\hat{r}}_y{\hat{r}}_z)}{{\mathrm{rV}}_S}\&CenterDot;{\stackrel{\&CenterDot;}{u}}_z(x^{\&prime;},y^{\&prime;},0,t+t_S)\end{array}\right]{\mathrm{dx}}^{\&prime;}{\mathrm{dy}}^{\&prime;}---\left(6\right)$

In the many component explorations of two dimension, use the harmomegathus focus usually, so only have P-P ripple and P-SV ripple on the theoretical seismogram.If making vertical is the z direction of principal axis, line direction is the x direction of principal axis, then calm field information on the y component.(5) formula, (6) formula are reduced to two-dimension elastic ripple Kirchhoff integral method vector shift formula.

$u_{\mathrm{PP}}(x,z,t)=\frac{1}{4\mathrm{\&pi;}}\underset{-\&infin;}{\overset{\&infin;}{\&Integral;}}[\frac{{2v}^2{\hat{r}}_x{\hat{r}}_z}{{\mathrm{rV}}_P}\&CenterDot;{\stackrel{\&CenterDot;}{u}}_x(x^{\&prime;},t+t_P)+\frac{1-{2v}^2{\hat{r}}_x{\hat{r}}_x}{{\mathrm{rV}}_P}\&CenterDot;{\stackrel{\&CenterDot;}{u}}_z(x^{\&prime;},t+t_P)]{\mathrm{dx}}^{\&prime;}---\left(7\right)$

$u_{\mathrm{PS}}(x,z,t)=\frac{1}{4\mathrm{\&pi;}}\underset{-\&infin;}{\overset{\&infin;}{\&Integral;}}[\frac{1-2{\hat{r}}_z{\hat{r}}_z}{{\mathrm{rV}}_S}\&CenterDot;{\stackrel{\&CenterDot;}{u}}_x(x^{\&prime;},t+t_S)+\frac{2{\hat{r}}_x{\hat{r}}_z}{{\mathrm{rV}}_S}\&CenterDot;{\stackrel{\&CenterDot;}{u}}_z(x^{\&prime;},t+t_S)]{\mathrm{dx}}^{\&prime;}---\left(8\right)$

Yet real medium is non-homogeneous and anisotropic, and underground diffraction is put the PP ripple of face of land acceptance point and the raypath of PS ripple no longer is a straight line, both raypaths, when walking and emergence angle all inequality.As shown in Figure 3, wherein dotted line is represented S wave ray path, and solid line is represented P wave ray path.Therefore need improve the homogeneous isotropic medium offset equation.T, t when walking _p, t _s, path l _p, l _s, emergence angle θ _p, θ _sWith equivalent p-and s-wave velocity V _pV _sRay tracing and velocity analysis by heterogeneous anisotropic media are tried to achieve respectively.

Therefore (7) formula, (8) formula can be improved

$u_{\mathrm{PP}}(x,z,t)=\frac{1}{4\mathrm{\&pi;}}\underset{-\&infin;}{\overset{\&infin;}{\&Integral;}}[\frac{{1-2v}^2{\sin }^2{\mathrm{\&theta;}}_P}{{l_{P}V}_P}\&CenterDot;{\stackrel{\&CenterDot;}{u}}_z(x^{\&prime;},t+t_P)+\frac{{2v}^2\sin {\mathrm{\&theta;}}_P{\cos \mathrm{\&theta;}}_P}{l_P{V}_P}\&CenterDot;{\stackrel{\&CenterDot;}{u}}_x(x^{\&prime;},t+t_P)]{\mathrm{dx}}^{\&prime;}---\left(9\right)$

$u_{\mathrm{PS}}(x,z,t)=\frac{1}{4\mathrm{\&pi;}}\underset{-\&infin;}{\overset{\&infin;}{\&Integral;}}[\frac{1-2\mathrm{co}{s}^2{\mathrm{\&theta;}}_S}{{l_{S}V}_S}\&CenterDot;{\stackrel{\&CenterDot;}{u}}_x(x^{\&prime;},t+t_S)+\frac{2\cos {\mathrm{\&theta;}}_S\sin {\mathrm{\&theta;}}_S}{{l_{S}V}_S}\&CenterDot;{\stackrel{\&CenterDot;}{u}}_z(x^{\&prime;},t+t_S)]{\mathrm{dx}}^{\&prime;}---\left(10\right)$

More than (9) formula, two equations of (10) formula be exactly the elastic wave Kirchhoff integral method vector shift formula that is applicable to heterogeneous anisotropic media.

Present homogeneous isotropism Kirchhoff integral offset algorithm (11) formula of contrast based on ACOUSTIC WAVE EQUATION:

$u(x,y,t)=-\frac{1}{2\mathrm{\&pi;}}\underset{A}{\text{\&Integral;\&Integral;}}\frac{\cos \mathrm{\&theta;}}{{\mathrm{lv}}_{\mathrm{mig}}}\stackrel{\&CenterDot;}{u}(x^{\&prime;},y^{\&prime;},0,t+\frac{1}{v_{\min }}){\mathrm{dx}}^{\&prime;}{\mathrm{dy}}^{\&prime;}---\left(11\right)$

(9) the vector elasticity wave field Kirchhoff integral offset algorithm in the non-homogeneous anisotropic medium of formula, (10) formula has following advantage:

1) has the characteristic of wave field separation, separate PP, PS ripple when can accomplish to be offset;

2) migration result is exported PP wave profile and PS wave profile respectively, is beneficial to the AVO feature of analyzing PP ripple and PS ripple respectively;

3) use the Green function of elastic oscillation equation, improved guarantor's width of cloth of algorithm;

4) be applicable to heterogeneous anisotropic media.

According to the said method principle, prestack depth migration method performing step of the present invention is as follows:

Step 1, geological data is carried out migration velocity analysis;

Step 2, when calculating away by ray tracing, arc length, emergence angle, incident angle;

Step 3, carry out migration aperture and calculate;

Step 4, utilize Kirchhoff integral method vector shift formula to carry out pre-stack depth migration.

Wherein, migration velocity analysis adopts common imaging point velocity analysis method based on deep focus in the step 1; In the step 2, ray casting when adopting direct rays back tracking method or employing to walk based on the minimum of Fermat principle; Adopt in the step 4 and be applicable to that the elastic wave Kirchhoff integral method vector shift formula of heterogeneous anisotropic media carries out pre-stack depth migration.

The flow process of prestack depth migration method performing step of the present invention is with reference to shown in Figure 4.

Wherein the flow process of vector pre-stack depth migration algorithm is as follows:

A. (geological data generally is that unit organizes with the seismic trace, adopts the SEG-Y stored in file format to the SEG-Y data.The SEG-Y form is one of standard reference tape data layout that is proposed by SEG (Society of Exploration Geophysicists, geophysics's exploration geophysicists), and it is the most general form of petroleum prospecting industry geological data) all roads circulations;

A) one SEG-Y data of input, and read shot point geophone station coordinate (judge it is X component or Z component during input, when different components will adopt different walking and Green function) by trace header;

B) with coarse travel timetable interpolation be and the identical grid of skew output section;

C) for underground all imaging point circulations;

● judge whether in PP ripple migration aperture;

If ● in pore diameter range, then calculate the Green function of this point, the pairing energy PP component of imaging point that is added to when reading this PP ripple and walking;

● judge whether in PS ripple migration aperture;

If ● in pore diameter range, then calculate the Green function of this point, the pairing energy PS component of imaging point that is added to when reading this PS ripple and walking;

D) to the imaging point loop ends;

The loop ends of B.SEG-Y data volume writes file with migration result, and skew finishes.

Its total algorithm realization flow is as shown in table 1.

The flow table of table 1 vector pre-stack depth migration algorithm

In the implementation procedure of prestack depth migration method of the present invention, the selection of migration aperture and the acquisition of migration velocity are two major control factors that are related to offset effect and counting yield, and the present invention has carried out following improvement to it.

Migration aperture

Select suitable migration aperture can improve the image quality that Kirchhoff integral method vector shift formula carries out pre-stack depth migration.Migration aperture is crossed conference and is reduced signal to noise ratio (S/N ratio); Otherwise then can omit effective information.Therefore PP ripple, the PS ripple raypath difference in medium needs to adopt different apertures.

The selection in aperture is the factor that has the greatest impact for imaging results (protecting the width of cloth, accuracy).The skew of Fresnel band, Gaussian beam skew, the finite aperture skew that are proposed in the existing literature in fact all are the actual positions of seeking the boundary reflection energy by ray theory, open up the aperture around this point.In this process, need to repeat the computing of trigonometric function and inverse trigonometric function, need very big calculated amount.And in the migration process, need disposable single big gun data of reading in.Therefore in actual production, all be based on the horizontal interface hypothesis, the perforate footpath when running into the inclination angle, just simply enlarges the aperture near CMP (Common Middle Point, common midpoint).

The method principle that the present invention adopts is identical with said method, but improves on algorithm, uses center point coordinate X _mReplace the coordinate of shot point and acceptance point respectively with offset distance h, avoided the calculating of trigonometric function and inverse trigonometric function, improved operation efficiency, and gone for offset method based on input channel.

Migration velocity analysis

Migration velocity analysis among the present invention adopts the common imaging point velocity analysis method based on deep focus, its basis is RCA method (Liu Z.and Bleistein N.Migration velocity analysisi:Theory and an iterative algorithm, Geophysics, 1995,60:142-153), this method is applicable to any geophone offset, inclination angle and velocity variations.In being total to image gather, imaging depth z is the function of geophone offset h.When migration velocity equaled underground true velocity, imaging depth and geophone offset were irrelevant, carry out Velocity Updating as criterion, and more new formula is as follows:

$\mathrm{\&delta;\&lambda;}=-\frac{{\mathrm{\&Sigma;}}_{k=1}^K{\mathrm{\&Sigma;}}_{j=1}^m(g_j^{\left(k\right)}-\stackrel{\&OverBar;}{{\hat{g}}^{\left(k\right)}})(z_j^{\left(k\right)}-\stackrel{\&OverBar;}{{\hat{z}}^{\left(k\right)}})}{{\mathrm{\&Sigma;}}_{k=1}^K{\mathrm{\&Sigma;}}_{j=1}^m{(g_j^{\left(k\right)}-\stackrel{\&OverBar;}{{\hat{g}}^{\left(k\right)}})}^2}$

Wherein, δ λ is the Velocity Updating step-length, and new iteration speed is V+ δ λ; The K representative is imaging point altogether, and M represents geophone offset; Wherein,

$g(x,h)=\frac{t_s+t_r}{\cos {\mathrm{\&theta;}}_s+{\cos \mathrm{\&theta;}}_r}$

Try to achieve by ray tracing, t is the interlayer whilst on tour;

The imaging depth of representing different geophone offsets; Z is by picking up acquisition on the image gather altogether.

Be the realization flow of migration velocity analysis as shown in Figure 5, its step comprises:

Step a1, initial geologic model, and enter step a2;

Step a2, set up initial velocity model, and enter step a3;

Step a3, carry out the PSDM of based target layer position bottom boundary, and enter step a4;

Step a4, carry out migration velocity analysis, and enter step a5;

Step a5, judged whether to set up complete rate pattern, if for "Yes" then enter step 2, if for "No" then enter step a6;

Step a6, carry out the PSDM of based target layer position bottom boundary, and enter step a7;

Step a7, pick up the destination layer bottom boundary, and enter step a8;

Step a8, renewal speed model move to down one deck position with destination layer, and enter step a3.

Wherein, complete rate pattern is described in the step a5: from shallow complete in speed, density, the layer position buried depth seismogeology model as parameter to dark each layer position.

This method fast convergence rate is applicable to various rate patterns and geophone offset.

Fig. 6 A is depicted as employing based on the PP ripple design sketch after the time migration after stack method of CCP road collection is offset 2D3C geological data of grand celebration altogether; Fig. 6 B is depicted as employing based on the PS ripple design sketch after the time migration after stack method of CCP road collection is offset 2D3C geological data of grand celebration altogether; Fig. 7 A is depicted as the PP ripple design sketch that adopts after the migration processing that the inventive method carries out a 2D3C geological data of grand celebration; Fig. 7 B is depicted as the PS ripple design sketch that adopts after the migration processing that the inventive method carries out a 2D3C geological data of grand celebration.Contrasted as can be seen by Fig. 6 A, Fig. 6 B, Fig. 7 A, Fig. 7 B, prestack depth migration method of the present invention is compared with previous method, has better effect:

The first, near surface (0-500m) imaging effect is good.

The second, imaging resolution improves.Shown in Fig. 6 A and Fig. 6 B, adopt, can only see it being many continuous lineups on the time migration section based in the resulting PP ripple and the effect of PS ripple after cheap of the time migration after stack method of CCP road collection altogether.And after adopting prestack depth migration method of the present invention, the resolution of lineups improves: shallow-layer (1000-2000m) lineups shown in Fig. 7 A and Fig. 7 B are more clear.It can also be seen that from Fig. 7 A and Fig. 7 B lineups are not continuous, but have a lot of little fractures, but in Fig. 6 A and Fig. 6 B, do not have this effect; The lineups dominant frequency at 3000 meters, deep of Fig. 7 A and Fig. 7 B obviously improves, and can see the tomography of can't see on two time migration sections.

Three, the continuity of the lineups at 3000 meters improves on the transformed wave imaging section, and it is more clear to construct.

Be illustrated in figure 8 as the simulated data model of prestack depth migration method of the present invention.Defining this model parameter is: 1000m * 1000m, four layers; Speed (m/s) Vp=2500,3000,3500,4000,4500; Vs=1500,1800,2100,2400,2700; Density (g/cm ³)=2.250,2.250,2.250,2.250,2.250.Adopt the elastic wave finite difference just drilling many big guns record; Recording geometry: track pitch 5m, 5 big guns, every big gun 200 roads, initial shot point coordinate x=100m, shot interval dx=200m; Sampling interval 0.5ms, sampling number: 2401.Fig. 9 A is the geological data Z component record diagram of the data model of Fig. 8, and Fig. 9 B is the geological data X component record diagram of the data model of Fig. 8.The migration result that adopts prestack depth migration method of the present invention is shown in Figure 10 (compressional wave migrated depth section), Figure 11 (converted shear wave migrated depth section); It shown in Figure 12 (compressional wave migrated depth section) and Figure 13 (converted shear wave migrated depth section) the ACOUSTIC WAVE EQUATION migration result that adopts with inclination angle perforate footpath; Shown in Figure 14 (compressional wave migrated depth section) and Figure 15 (converted shear wave migrated depth section) result who adopts the skew of constant aperture ACOUSTIC WAVE EQUATION.Imaging point grid 100*100, dx=10m, dz=10m; Migration aperture aper=0.4, promptly aperture radius takes advantage of 0.4 for the imaging point depth z.

Contrast as can be seen from Figure 10-Figure 15, the offset effect of the inventive method is best, and the structure playback is accurate, clear; Secondly be ACOUSTIC WAVE EQUATION offset method by inclination angle perforate footpath; The poorest is constant aperture ACOUSTIC WAVE EQUATION offset method.

The above only is preferred embodiment of the present invention, and is in order to restriction the present invention, within the spirit and principles in the present invention not all, any modification of being done, is equal to replacement, improvement etc., all should be included within protection scope of the present invention.

Claims (5)

1. prestack depth migration method may further comprise the steps:

Step 1, geological data is carried out migration velocity analysis;

Step 2, when calculating away by ray tracing, arc length, emergence angle, incident angle;

Step 3, carry out migration aperture and calculate;

Step 4, utilize Kirchhoff integral method vector shift formula to carry out pre-stack depth migration.

2. prestack depth migration method according to claim 1 is characterized in that: migration velocity analysis described in the step 1 adopts the common imaging point velocity analysis method based on deep focus.

3. prestack depth migration method according to claim 2 is characterized in that described step 1 comprises the steps:

Step a1, initial geologic model, and enter step a2;

Step a2, set up initial velocity model, and enter step a3;

Step a3, carry out the pre-stack depth migration of based target layer position bottom boundary, and enter step a4;

Step a4, carry out migration velocity analysis, and enter step a5;

Step a5, judged whether to set up complete rate pattern, if for "Yes" then enter step 2, if for "No" then enter step a6;

Step a6, carry out the pre-stack depth migration of based target layer position bottom boundary, and enter step a7;

Step a7, pick up the destination layer bottom boundary, and enter step a8;

Step a8, renewal speed model move to down one deck position with destination layer, and enter step a3.

4. prestack depth migration method according to claim 1 is characterized in that: in the step 2, and ray casting when adopting direct rays back tracking method or employing to walk based on the minimum of Fermat principle.

5. prestack depth migration method according to claim 1 is characterized in that: adopt in the step 4 to be applicable to that the elastic wave Kirchhoff integral method vector shift formula of heterogeneous anisotropic media carries out pre-stack depth migration.

Images