CN104101901A - Converted-wave curved ray amplitude-reserved anisotropic pre-stack time offset time method - Google Patents

Abstract

The invention provides a converted-wave curved ray amplitude-reserved anisotropic pre-stack time offset time method, belonging to the geophysical exploration field. The method comprises steps of inputting three-dimensional converted-wave pre-stack gather record and a converted-wave anisotropic velocity body, beginning imaging, calculating longitudinal wave offset speed and transverse wave speed offset speed at the imaging sample point in series, calculating longitudinal wave anisotropic parameters and transverse wave anisotropic parameters, calculating longitudinal wave travel time and transverse wave travel time, calculating converted wave anisotropic travel time, calculating a cosine of a ray angle of incidence, calculating a cosine of a ray angle of emergence, calculating geometric diffusion term, calculating amplitude weighting coefficient, finishing amplitude preserving migration processing, and performing imaging by performing the above processing on all imaging points.

Description

A kind of transformed wave curved rays is protected width anisotropy prestack time migration method

Technical field

The invention belongs to geophysical survey field, be specifically related to a kind of transformed wave curved rays and protect width anisotropy prestack time migration method.

Background technology

While producing artificial earthquake wave field on earth's surface with vertical epicenter excitation, wherein a part of compressional wave energy is propagated downwards, arrive subsurface reflective boundary, on the interface in underground wave impedance reflection horizon, incident longitudinal wave is except occurring reflected P-wave (PP), according to Zuo Pulizi equation, can also be converted to shear wave (S), direction is propagated earthward.Receive with three-component seismometer record on earth's surface, just can obtain compressional wave (PP) and transformed wave (PS) multi-component seismic data.PS ripple pre-stack time migration (pre-stack time migration) formation method, can solve the underground structure imaging problem of special reservoir, as serious in PP ripple attenuation by absorption under gas cloud, cause reservoir structure imaging fuzzy, and PS ripple is less because S speed is affected by fluid generally, depend on rock skeleton speed, irrelevant with the contained fluid of rock skeleton hole.Therefore can obtain reservoir structure imaging under gas cloud.In PS ripple data collection, in the switching energy of offset distance far away stronger, therefore in multi-component seismic data acquisition design the observation program of the normal large offseting distance adopting.Due to the lamellarity of medium, a kind of anisotropy (V TI) properties influence with vertical axis of symmetry each level same sex medium compressional wave and shear transit time simultaneously, because p-and s-wave velocity has directivity in the time that this medium is propagated.The reflection coefficient of seismic event and wave trajectory and ray angle of departure are approximation relation.Based on the migration algorithm of kirchhoff iterated integral type, need a rational amplitude weight function to carry out correction of amplitude to the true amplitude recovery of seismic event, keep relative real amplitude information.

At present in the whilst on tour computing method of transformed wave anisotropy (VTI) medium, the most classical basic representation is that Thomsen L. in 1999 is published in the U.S. " geophysics " " Converted-wave reflection seismology over inhomogeneous, anisotropic media " literary composition, the method in the time that downgoing P-wave is walked, calculate in except velocity of longitudinal wave parameter, also has a compressional wave anisotropic parameters, in calculating in the time that up S ripple is walked, except shear wave velocity, also consider shear wave anisotropic parameters, and at application United States Patent (USP) in 2000.Citation form still adds anisotropic parameters in hyperbolic curve whilst on tour formula.

About non-double curve line style seimic travel time computing formula, can trace back to the earliest " Velocity spetradigital computer derivation and applications of velocity function " literary composition that Taner M.T. in 1969 is published in the U.S. " geophysics ", elaborate in this article the hyperbolic curve Computation of seismic travel time formula that replaces individual layer root-mean-square velocity by layering and interval velocity computing formula.Constantly there is the seismic event T-X curve computing formula of new non-double curve line style to release thereupon, " the A new travel time estimation method for horizontal strata " paper and the E.Blias in 2007 that deliver in (SEG) annual meeting of U.S. geophysical survey person association as nearest 2005 M.Turhan and Taner learn Canadian geophysical survey person " Long-Offset Approximations to NMO Function for a Layered Subsurface " paper that (CSEG) annual meeting is delivered, and are all to improve constantly Computation of seismic travel time precision and counting yield.

The most representative method of transformed wave anisotropy prestack time migration method is to be published in anisotropy project group (EAP) Li Xiangyang of geologic examination office of Britain and Dai Hengchang calendar year 2001 " Anisotropy migration and model building for 4C seismic data " paper in U.S. SEG meeting, this article has been set forth the processing thinking that PS ripple is equivalent to C ripple, introduce the method for C ripple anisotropy velocity modeling, as vertical speed ratio, effective velocity ratio, the computing method of C ripple root-mean-square velocity and C ripple anisotropic parameters, in PS ripple prestack time migration method, how to convert C ripple anisotropy rate pattern to P ripple and S wave propagation velocity and anisotropic model, thereby realize PS ripple anisotropy pre-stack time migration disposal route.In the method, when whilst on tour, computing formula has still adopted hyperbolic curve form, has just added anisotropic parameters.And amplitude weight function calculates the simple time, the relevant weight coefficient of speed of adopting, irrelevant with the geometric shape of x-ray angle and recording geometry.

SUN and the Martinez etc. of PGS company of the U.S., in U.S. SEG in 2003 meeting, deliver " 3DKirchhhoff PS-wave prestack time migration for V (z) and VTI media " paper, the high-order that provides first offset distance 6 powers in whilst on tour calculates calculates whilst on tour formula, and on amplitude weight function, has considered angle of departure, the recording geometry geometric relationship of ray.And with within 2004, applied for United States Patent (USP).The LEE Seongbok of the said firm in 2006 has applied for " prestack migration method of amplitude reserving " international monopoly technology, and patent has disclosed the method for the amplitude reserving of kirchhoff time migration.This algorithm has comprised by the constant speed estimation geometrical attenuation factor, according to rate pattern accurate Calculation incident angle and the emergence angle of time domain.Estimation weighting factor and angle of departure carry out the calculating of imaging point.

The PS ripple prestack time migration method of current domestic employing, substantially be by technology introduction, the prestack time migration method of anisotropy project group (EAP) Li Xiangyang of geologic examination office of Britain and Dai Hengchang is transplanted to a processing module in the disposal system that becomes independent research, but it is the earthquake travel-time equation based on hyperbolic-type that its whilst on tour calculates, only add VTI anisotropic parameters at the earthquake travel-time equation of hyperbolic-type, amplitude weight function does not have yet can Consideration of Three-dimensional geometric shape and the influence factor such as ray take-off angle, be referred to as the anisotropy prestack time migration method of the non-guarantor's width of straight line.

At present in C ripple prestack time migration method, about the computing formula of transformed wave whilst on tour be:

$t_c=\sqrt{{\left(\frac{t_{c0}}{1+{\mathrm{\γ}}_0}\right)}^2+\frac{{(x+h)}^2}{V_{\mathrm{pn}}^2}-{2\mathrm{\η}}_{\mathrm{eff}}{\mathrm{\Δt}}_P^2}+\sqrt{{\left(\frac{{\mathrm{\γ}}_0{t}_{c0}}{1+{\mathrm{\γ}}_0}\right)}^2+\frac{{(x-h)}^2}{V_{\mathrm{sn}}^2}+2{\mathrm{\ζ}}_{\mathrm{eff}}{\mathrm{\Δt}}_s^2}---\left(1\right)$

T in formula _cfor C ripple whilst on tour, γ ₀for compressional wave and shear wave vertical speed ratio, t _c0for C ripple t ₀time, x is offset distance, h is imaging point horizontal level, V _pnand V _snbe respectively compressional wave and shear wave stack velocity, η _effand ξ _effbe respectively compressional wave and shear wave anisotropic parameters.Δ t _pwith Δ t _sdetermined by following formula:

${\mathrm{\Δt}}_P^2=\frac{{(x+h)}^4}{V_{\mathrm{pn}}^2[t_{p0}^2{V}_{\mathrm{pn}}^2+(1-{2\mathrm{\η}}_{\mathrm{eff}}){(x+h)}^2]}$

${\mathrm{\Δt}}_s^2=\frac{{(x-h)}^4}{V_{\mathrm{sn}}^2[t_{s0}^2{V}_{\mathrm{sn}}^2+{(x-h)}^2]}$

(2)

The form of weighting function is:

$w=\frac{1}{\sqrt{t_c{V}_{\mathrm{pn}}{V}_{\mathrm{sn}}}}---\left(3\right)$

In the method, when whilst on tour, computing formula has still adopted hyperbolic curve form, has just added anisotropic parameters.And amplitude weight function calculates simple time of employing, the relevant weight coefficient of speed.

But, apply this hyperbolic-type whilst on tour computing formula (1), carry out the processing of kirchhoff pre-stack time migration, can cause the skew of large offseting distance data to occur inclined to one side phenomenon (seeing Fig. 1 square frame).Carry out amplitude of deflection correction with this simple distance weighting function formula (3), also weighting phenomenon was appearred in the road collection that causes large offseting distance, be equal to without amplitude weight and proofread and correct (seeing Fig. 2 square frame)

Summary of the invention

The object of the invention is to solve the difficult problem existing in above-mentioned prior art, provide a kind of transformed wave curved rays to protect width anisotropy prestack time migration method, can improve converted-wave prestack time migration imaging reflection wave groups energy focusing, overcome low frequency phenomenon; Improve migration precision, make breakpoint clear; Enhanced Imaging section signal to noise ratio (S/N ratio); Imaging road collection (CIP) of output meets amplitude variation with Offset (AVO) rule more.Final migrated section contributes to structure elucidation, output CIP road collection to be applicable to AVO analysis and AVO inverting application

The present invention is achieved by the following technical solutions:

Said method comprising the steps of:

(1) input three-dimensional transformed wave prestack road collection record and transformed wave anisotropy body of velocity, comprise that ripple vertical speed compares γ in length and breadth ₀, velocity equivalent compares γ _eff, transformed wave anisotropic parameters x _effwith transformed wave migration velocity V _cnthese four parameters;

(2) start a track data of input to carry out imaging processing;

(3) four parameters of applying step (1) input, are calculated to be the compressional wave migration velocity V of decent some place (sampling point circulates successively, finishes sampling point from calculating initial sampling point Dao Zhe road) successively _pnwith shear wave migration velocity V _sn;

Calculate compressional wave anisotropic parameters η _effwith shear wave anisotropic parameters ξ _eff;

Calculate compressional wave whilst on tour t _p6with shear wave whilst on tour t _s6;

(4) the compressional wave curved rays whilst on tour t that applying step (3) is tried to achieve _p6, shear wave whilst on tour t _s6, compressional wave anisotropic parameters η _effwith shear wave anisotropic parameters ξ _effcalculate transformed wave anisotropy whilst on tour t _c;

(5) according to compressional wave migration velocity and compressional wave curved rays whilst on tour t _p6calculate ray incident angle α _scosine cos (α _s); According to shear wave migration velocity and shear wave curved rays whilst on tour t _s6calculate ray emergence angle α _gcosine cos (α _g);

(6) computational geometry diffusion term, L and sin (α);

(7) according to transformed wave amplitude weight function, calculated amplitude weighting coefficient W (ξ, R);

(8) according to transformed wave offset equation, application amplitude weight coefficient W (ξ, R) and transformed wave whilst on tour t _cand transformed wave input wave field derivative term, complete transformed wave curved rays and protect width migration processing (just referring to step [12] below);

(9) determine whether last imaging point, if so, proceed to step (10), if not, return to step (2);

(10) one imaging finish.

In described step (3), be the compressional wave migration velocity V that is calculated to be respectively picture point place according to formula (4) _pnwith shear wave migration velocity V _sn:

$V_{\mathrm{pn}}^2=V_{\mathrm{cn}}^2\frac{{\mathrm{\γ}}_{\mathrm{eff}}(1+{\mathrm{\γ}}_0)}{1+{\mathrm{\γ}}_{\mathrm{eff}}}$

(4)

$V_{\mathrm{sn}}^2=V_{\mathrm{cn}}^2\frac{(1+{\mathrm{\γ}}_0)}{{\mathrm{\γ}}_0(1+{\mathrm{\γ}}_{\mathrm{eff}})}$

In described step (3), be to be calculated to be respectively the compressional wave anisotropic parameters η of picture point place according to formula (5) _effwith shear wave anisotropic parameters ξ _eff:

${\mathrm{\η}}_{\mathrm{eff}}=\frac{{\mathrm{\χ}}_{\mathrm{eff}}}{({\mathrm{\γ}}_0-1){\mathrm{\γ}}_{\mathrm{eff}}^2}$

(5)

${\mathrm{\ζ}}_{\mathrm{eff}}=\frac{{\mathrm{\χ}}_{\mathrm{eff}}}{{\mathrm{\γ}}_0-1}$

In described step (3), be to calculate compressional wave whilst on tour t according to formula (10) and (11) _p6with shear wave whilst on tour t _s6:

$t_{p6}=t_{p4}(1+\frac{1}{2}k\frac{c_{4p}{x}_p^6}{t_{p4}^2})---\left(10\right)$

$t_{s6}=t_{s4}(1+\frac{1}{2}k\frac{c_{4s}{x}_s^6}{t_{s4}^2})---\left(11\right)$

Wherein, k is high-order small quantity, as constant coefficient item, and value 0.0～2.0;

T _p4for imaging point place compressional wave 4 rank item whilst on tours, t _s4for imaging point place shear wave 4 rank item whilst on tours;

$t_{p4}=\sqrt{c_{1p}+c_{2p}{x}_p^2+c_{3p}{x}_p^4}---\left(8\right)$

$t_{s4}=\sqrt{c_{1s}+c_{2s}{x}_s^2+c_{3s}{x}_s^4}---\left(9\right)$

Wherein, x _prepresent the distance from focal point to imaging point, x _srepresent the distance from geophone station to imaging point;

$c_1=a_1^2=\frac{2z}{v_r},$ $c_2=\frac{a_1}{a_2}=\frac{1}{{v_r}^2},$ $c_3=\frac{a_2^2-a_1{a}_3}{4a_2^4},$ $c_4=\frac{{2a}_1{a}_3^2-a_1{a}_2{a}_4-a_2^3{a}_3}{{8a}_2^7}---\left(6\right)$

$a_i=2\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{v}_i^{2i-3}{h}_k---\left(7\right)$

V in formula _i, h _kbe respectively interval velocity and zone thickness, n is hierarchy number.

C _1p, c _2p, c _3pand c _4pfor the coefficient of compressional wave, adopt P-wave interval velocity to utilize formula (6) to calculate; c _1s, c _2s, c _3sand c _4sfor the coefficient of shear wave, adopt S-wave interval velocity to utilize formula (6) to calculate.

Described step (4) is to calculate transformed wave anisotropy whilst on tour t according to formula (12) _c:

$t_c=\sqrt{t_{p6}^2-2{\mathrm{\η}}_{\mathrm{eff}}{\mathrm{\Δt}}_p^2}+\sqrt{t_{s6}^2+{2\mathrm{\ζ}}_{\mathrm{eff}}{\mathrm{\Δt}}_s^2}---\left(12\right)$

Δ t in formula _pwith Δ t _scalculated by formula (2):

$\mathrm{\Δ}{t}_P^2=\frac{{(x+h)}^4}{V_{\mathrm{pn}}^2[t_{p0}^2{V}_{\mathrm{pn}}^2+(1+2{\mathrm{\η}}_{\mathrm{eff}}){(x+h)}^2]}$

(2)。

$\mathrm{\Δ}{t}_s^2=\frac{{(x-h)}^4}{V_{\mathrm{sn}}^2[t_{s0}^2{V}_{\mathrm{sn}}^2+{(x-h)}^2]}$

Described step (5) utilizes formula (13) to calculate:

$\cos \left({\mathrm{\α}}_s\right)=\sqrt{1-{\left(\frac{{\∂T}_p}{\∂x_p}\right)}^2{V_p}^2}$

$\cos \left({\mathrm{\α}}_g\right)=\sqrt{1-{\left(\frac{{\∂T}_s}{\∂x_s}\right)}^2{V_s}^2}---\left(13\right).$

Tp in above formula, that Ts adopts is compressional wave whilst on tour t _p6with shear wave whilst on tour t _s6.

Described step (6) is achieved in that

First calculate offset midpoint position M (Xm _,y _m) coordinate:

$X_m=\frac{X_s+X_g}{2}$

(14)

$Y_m=\frac{Y_s+Y_g}{2}$

X in formula _sand X _gbe respectively big gun inspection coordinate;

Then calculate reflection spot at ground projected position Q (X _q, Y _q) coordinate:

X _q＝CDP _x-X _m

Y _q＝CDP _y-Y _m (15)

CDP in formula _x, CDP _yfor the road header of present input data; L is the length of QM, is shown with end points coordinates table:

$L=\sqrt{{\left(X_q\right)}^2+{\left(Y_q\right)}^2}---\left(16\right)$

And the sine of the angle α of two straight line MQ and SG can calculate with the slope of two straight lines:

$\sin \mathrm{\α}=\frac{k_2-k_1}{\sqrt{(1+{k_1}^2)(1+{k_2}^2)}}---\left(17\right)$

Described step (7) is to adopt formula (18) to calculate:

$W(\mathrm{\ξ},R)=\frac{\sqrt{\cos {\mathrm{\α}}_s\cos {\mathrm{\α}}_g}}{V_{p0}}\left\{T_r\right[\frac{T_s}{T_p^2}+\frac{T_p}{T_s^2}]+\frac{4L^2{H}^2{\sin }^2\left(\mathrm{\α}\right)}{T_p{T}_s{T}_r^2{V}_c^4}\}---\left(18\right)$

T in formula _rfor the time at imaging point place, T _pand T _s, the compressional wave whilst on tour t that adopts step (3) to calculate _p6with shear wave whilst on tour t _s6; V _cfor the migration velocity of imaging point, adopt described transformed wave migration velocity V _cn,, V _p0be compressional wave near-surface velocity, H is geophone offset.

Compared with prior art, the invention has the beneficial effects as follows: the high-order anisotropy whilst on tour computing formula of application the inventive method, on skew output CIP road collection, showing centering offset gather far away has obvious offset correction to even up effect.Therefore relatively existing method, the inventive method can be improved converted-wave prestack time migration imaging reflection wave groups energy focusing, overcomes low frequency phenomenon, and tracing of horizons is explained reliable; Improve migration precision, make breakpoint clear; The weighting function of application this method, Enhanced Imaging section signal to noise ratio (S/N ratio), output CIP road collection meets amplitude variation with Offset (AVO) rule more.Improve generally imaging section quality, contribute to profile construction to explain, output CIP road collection is applicable to AVO analysis and AVO inverting application.

Aspect seismic structure imaging and seismic properties inverting, apply this technology and can complete transformed wave three peacekeeping two dimension pre-stack time migrations.Result can improve pre-stack time migration section image quality, strengthens the signal to noise ratio (S/N ratio) of migrated section, and deflection focusing is good, and tracing of horizons is explained reliable, and migration is accurate, and breakpoint is clear, and substratum resolution is high, and the imaging road collection CIP of output has AVO feature.Migrated section contributes to structure elucidation, CIP road collection to be applicable to prestack AVO inverting, therefore in reservoir exploration and exploitation, has a good application prospect.

Brief description of the drawings

Fig. 1 be imaging road collection (CIP) process a long way occurred inclined to one side.

Tu2Shi CIP road collection is processed a long way and was occurred weighting.

Fig. 3 (a) is Depth Domain rate pattern.

Fig. 3 (b) is Wave equation forward modeling big gun record

Layered medium curved rays path, Fig. 4-1 schematic diagram.

Single-layer medium direct rays path, Fig. 4-2 schematic diagram.

Parameter-definition in Fig. 5 formula (12), R imaging point position, S focal point, G geophone station, Q be imaging point R in ground projection, M is that big gun is in cautious.H is geophone offset.

The three-dimensional PS ripple anisotropy pre-stack time migration CIP of the existing method direct projection in Fig. 6-1 collimation method road collection.

The CIP of the three-dimensional PS ripple of Fig. 6-2 patented method curved rays of the present invention method anisotropy pre-stack time migration.

The three-dimensional PS ripple of the existing method direct projection of Fig. 7 (a) collimation method anisotropy pre-stack time migration migrated section.

The three-dimensional PS ripple of Fig. 7 (b) patented method curved rays of the present invention method anisotropy pre-stack time migration migrated section.

Fig. 8 (a) processes CIP road collection without weighting function.

Fig. 8 (b) is that the simple amplitude weight function of existing method is processed CIP road collection.

Fig. 8 (c) is patented method weighting function formula of the present invention (12) processing CIP road collection.

The pre-stack time migration section of the non-guarantor's width of the existing method direct rays in Fig. 9-1.

Fig. 9-2 the inventive method curved rays is protected the pre-stack time migration processing profiles of width.

Figure 10 is the step block diagram of the inventive method.

Embodiment

Below in conjunction with accompanying drawing, the present invention is described in further detail: below in conjunction with accompanying drawing, the present invention is described in further detail:

Transformed wave curved rays of the present invention is protected width anisotropy prestack time migration method and is adopted first higher order term (offset distance 6 powers) anisotropy PS ripple non-double curve line style whilst on tour computing formula at home and abroad; in amplitude weight function calculates, having considered to gather geometric relationship and ray incident angle and the reflection angle calculating of recording geometry, is the transformed wave anisotropy prestack time migration method of the kirchhoff integral algorithm of a kind of high precision whilst on tour calculating and amplitude reserving.Computing formula of the present invention combines the advantage of EAP project team converted-wave prestack time migration method and PGS company converted-wave prestack time migration patented method, forms own unique algorithm.This method has arrived current international most advanced level.

The present invention is in the Computation of seismic travel time formula of Kirchhoff migration; adopt 6 rank items and vertical axis of symmetry isotropy (VTI) parameter of offset distance; improve transformed wave whilst on tour computational accuracy; more approach the actual seismic wave propagation time, the migration processing of the multi-component seismic data that in adaptation, offset distance far away gathers.In amplitude of deflection weighting function calculates, according to the rate pattern estimation three-dimensional geometry invasin of time domain, the incident angle of accurate Calculation ray and emergence angle, estimation weighting factor carries out correction of amplitude to each imaging point, recovers seismic event real amplitude information.

As shown in figure 10, the inventive method comprises the following steps:

[1] the pretreated three-dimensional transformed wave prestack of prestack road collection data were done in input, the data of inputting be must through static correction and residual static correction processing, prestack denoising particularly maximum value pick processing, data can be ACP (asymptotic line altogether) or CCP (transfer point altogether) road collection data, also can be shot gather data, input data possess normal road header (surveying wire size, the cautious coordinate of big gun, offset distance etc.).Input the transformed wave anisotropic three-dimensional body of velocity of transformed wave anisotropy migration velocity analysis modeling simultaneously, ripple vertical speed ratio, effective velocity ratio, transformed wave anisotropic parameters and 4 parameters of transformed wave migration velocity are in length and breadth provided respectively.

[2] the three-dimensional transformed wave anisotropy body of velocity of applying step [1] input, i.e. transformed wave migration velocity V _cn, vertical speed compares γ ₀compare γ with effective velocity _eff, be calculated to be respectively the compressional wave migration velocity V of picture point place according to formula (4) _pnwith shear wave migration velocity V _sn.

$V_{\mathrm{pn}}^2=V_{\mathrm{cn}}^2\frac{{\mathrm{\γ}}_{\mathrm{eff}}(1+{\mathrm{\γ}}_0)}{1+{\mathrm{\γ}}_{\mathrm{eff}}}$

(4)

$V_{\mathrm{sn}}^2=V_{\mathrm{cn}}^2\frac{(1+{\mathrm{\γ}}_0)}{{\mathrm{\γ}}_0(1+{\mathrm{\γ}}_{\mathrm{eff}})}$

[3] the three-dimensional transformed wave anisotropy body of velocity of applying step [1], i.e. transformed wave anisotropic parameters x _eff, vertical speed compares γ ₀compare γ with effective velocity _eff, be calculated to be respectively the compressional wave anisotropic parameters η of picture point place according to formula (5) _effwith shear wave anisotropic parameters ξ _eff.

${\mathrm{\η}}_{\mathrm{eff}}=\frac{{\mathrm{\χ}}_{\mathrm{eff}}}{({\mathrm{\γ}}_0-1){\mathrm{\γ}}_{\mathrm{eff}}^2}$

(5)

${\mathrm{\ζ}}_{\mathrm{eff}}=\frac{{\mathrm{\χ}}_{\mathrm{eff}}}{{\mathrm{\γ}}_0-1}$

The elementary tactics of the curved rays [4] showing according to Fig. 4-1 and Fig. 4-2 and direct rays layering and interval velocity, direct rays can equivalence be regarded one deck (Fig. 4-2) as from earth's surface to imaging point position, and speed is got root-mean-square velocity V _rms.Curved rays is divided into some layers of h by earth's surface to imaging point position _k(Fig. 4-1), the speed of every one deck is got speed V in its layer _i.Compressional wave and shear wave root-mean-square velocity field that applying step [2] calculates, estimate compressional wave and the S-wave interval velocity of each layer, how much consuming time relevant with calculating with computational accuracy (approximate true whilst on tour) layering is, and layering is more, computational accuracy is higher, and the time of calculating consumption is just more.It should be noted that: the initial time that requires input transformed wave rate pattern is zero or close to zero, otherwise affects computational accuracy.

[5] without loss of generality, high-order curved rays coefficient formulas is (6) and (7) formula, the lift height, layering number and the vertical S-wave interval velocity that utilize step [4] to obtain, respectively design factor c ₁, c ₂, c ₃and c ₄, in the time adopting P-wave interval velocity, we obtain the c of compressional wave _1p, c _2p, c _3pand c _4pcoefficient, in the time adopting S-wave interval velocity, obtains the c of shear wave _1s, c _2s, c _3sand c _4scoefficient,

$c_1=a_1^2=\frac{2z}{v_r},$ $c_2=\frac{a_1}{a_2}=\frac{1}{{v_r}^2},$ $c_3=\frac{a_2^2-a_1{a}_3}{{4a}_2^4},$ $c_4=\frac{{2a}_1{a}_3^2-a_1{a}_2{a}_4-a_2^3-a_3}{{8a}_2^7}---\left(6\right)$

In formula

$a_i=2\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{v}_i^{2i-3}{h}_k---\left(7\right)$

V in formula _i, h _kbe respectively interval velocity and zone thickness, n is hierarchy number.

[6] compressional wave c step [5] being calculated _1p, c _2p, c _3pand c _4pthe c of coefficient and shear wave _1s, c _2s, c _3sand c _4scoefficient, according to curved rays whilst on tour formula (8) and (9), is calculated to be respectively picture point place compressional wave 4 rank item whilst on tour t _p4, shear wave 4 rank item whilst on tour t _s4.

$t_{p4}=\sqrt{c_{1p}+c_{2p}{x}_p^2+c_{3p}+x_p^4}---\left(8\right)$

$t_{s4}=\sqrt{c_{1s}+c_{2s}{x}_s^2+c_{3s}{x}_s^4}---\left(9\right)$

X in formula _prepresent the distance from focal point to imaging point, x _srepresent the distance from geophone station to imaging point.According to curved rays whilst on tour formula (10) and (11), be calculated to be respectively picture point place compressional wave 6 rank item whilst on tour t _p6, shear wave 6 rank item whilst on tour t _s6.

$t_{p6}=t_{p4}(1+\frac{1}{2}k\frac{c_{4p}{x}_p^6}{t_{p4}^2})---\left(10\right)$

$t_{s6}=t_{s4}(1+\frac{1}{2}k\frac{c_{4s}{x}_s^6}{t_{s4}^2})---\left(11\right)$

In formula, k is high-order small quantity, as constant coefficient item, and value 0.0～2.0.Thereby complete at the locational descending compressional wave whilst on tour t of imaging point _p6with up shear wave whilst on tour t _s6calculate.

[7] calculate transformed wave curved rays anisotropy whilst on tour, the curved rays whilst on tour result that applying step [6] calculates, and step [3] compressional wave and the shear wave anisotropic parameters that calculate respectively, according to formula (12), (the anisotropy whilst on tour not refering in particular to is general term to calculate the anisotropy transformed wave whilst on tour of the curved rays on image space, the anisotropy whilst on tour that adds curved rays refers in particular to (12) formula, because use higher order term) t _c.

$t_c=\sqrt{t_{p6}^2-2{\mathrm{\η}}_{\mathrm{eff}}\mathrm{\Δ}{t}_p^2}+\sqrt{t_{s6}^2+2{\mathrm{\ζ}}_{\mathrm{eff}}\mathrm{\Δ}{t}_s^2}---\left(12\right)$

Δ t in formula _pwith Δ t _scalculate gained by formula (2).

Fig. 3 (a) and Fig. 3 (b) are direct rays and curved rays whilst on tour computational accuracy comparison diagram.Fig. 3 (a) is interval velocity model, and model is divided into 9 layers, and maximum interface depth reaches 4600 meters, interval velocity distribution range 1500-4800 meter per second.Utilize Wave equation forward modeling list big gun compressional wave record to be presented at Fig. 3 (b), the reflection wave hyperbolic curve of certain one deck in Fig. 3 (b) (1800 milliseconds of left and right) as figure in canescence show.The T-X curve that when two kinds of simultaneously superimposed demonstrations are walked, computing method obtain.Dark line is direct rays travel time curve, and light grey line is curved rays travel time curve.On travel time curve a long way, grayish curved rays travel time curve shows and more approaches linen wave equation travel time curve, and computational accuracy obviously improves, and the dark line of direct rays and wave equation travel time curve error obviously increase.

[8] protect width migration processing in order to realize transformed wave, the present invention has adopted the form of the weighting function of geometric relationship, accurate Calculation ray angle of departure and the source wavelet of estimation recording geometry.Formula (10) and (11) are asked about x _pand x _spartial derivative the p-and s-wave velocity field that obtains of applying step [2], calculate respectively the incident angle α of ray according to formula (13) _s(shot point place) and emergence angle α _gthe cosine form at (geophone station place).

$\cos \left({\mathrm{\α}}_s\right)=\sqrt{1-{\left(\frac{\∂T_p}{\∂x_p}\right)}^2{V_p}^2}$

$\cos \left({\mathrm{\α}}_g\right)=\sqrt{1-{\left(\frac{\∂T_s}{\∂x_s}\right)}^2}{V_s}^2---\left(13\right)$

[9] geometrical attenuation item calculates L and sin (α) calculating, first calculates offset midpoint position M (X _m, Y _m) coordinate (and referring to Fig. 5, R is imaging point position, and S is focal point, and G is geophone station, Q be imaging point R in ground projection, M is that big gun is in cautious.H is geophone offset.)：

$X_m=\frac{X_s+X_g}{2}$

(14)

$Y_m=\frac{Y_s+Y_g}{2}$

X in formula _sand X _gbe respectively big gun inspection coordinate (from seismic trace header).Then calculate reflection spot at ground projected position Q (X _q, Y _q) coordinate:

X _q＝CDP _x-X _m

(15)

Y _q＝CDP _y-Y _m

CDP in formula _x, CDP _yfor the road header of present input data.L is the length of QM, is shown with end points coordinates table:

$L=\sqrt{{\left(X_q\right)}^2+{\left(Y_q\right)}^2}---\left(16\right)$

And the sine of the angle α of two straight line MQ and SG can calculate with the slope of two straight lines

$\sin \mathrm{\α}=\frac{k_2-k_1}{\sqrt{(1+{k_1}^2)(1+{k_2}^2)}}---\left(17\right)$

[10] whole weighting function numerical procedure adopts the expression formula of formula (18),

$W(\mathrm{\ξ},R)=\frac{\sqrt{\cos {\mathrm{\α}}_s\cos {\mathrm{\α}}_g}}{V_{p0}}\left\{T_r\right[\frac{T_s}{T_p^2}+\frac{T_p}{T_s^2}]+\frac{{4L}^2{H}^2{\sin }^2\left(\mathrm{\α}\right)}{T_p{T}_s{T}_r^2{V}_c^4}\}---\left(18\right)$

T in formula _rfor the time at imaging point place, T _pand T _s, be the whilst on tour that step [6] is calculated, the whilst on tour time (shear wave time) that is respectively imaging point to whilst on tour time (compressional wave time) of shot point and imaging point to geophone station.V _cfor the migration velocity (being the transformed wave velocity field of step [1] input) of imaging point, C _p0be compressional wave near-surface velocity (the starting velocity field that the compressional wave time of calculating for step [2] is zero), H is geophone offset (from input channel header), and L and sin (α) have been calculated by step [9].Formula (18) has been described the impact of source wavelet, reflection and geometrical attenuation successively.

[11] application of weighting coefficient: the use of real amplitude weighting coefficient, need to be in migration aperture (given one group of time, angle pair: as 0.0:20.0; 3.0:40; 6.0:60) in limit, in the time that oval amplitude value (being the total whilst on tour of time/transformed wave of image space) is less than migration aperture, the amplitude weight coefficient that step [10] is calculated is set to zero, in the time that oval amplitude value is greater than migration aperture, weighting coefficient is multiplied by actual wave field amplitude, realizes amplitude compensation.

[12] Kirchhoff Summation Method of Migration method: the kirchhoff algorithm of pre-stack time migration is the process that the amplitude of seismic wave field is weighted to summation along diffraction curve, can be write as with integral expression:

$I(\mathrm{\τ},y,h)=\∫W(\mathrm{\τ},y,b,h)\frac{\∂}{\∂t}u(\mathrm{\τ},y,b,h)\mathrm{db}---\left(19\right)$

In formula, h is 1/2nd geophone offset, y is common midpoint coordinate, b is the distance (by trace header information acquisition) that Diffraction Imaging point departs from central point, and W is weighting function (step [10] calculates), and I is that transformed wave whilst on tour τ (is t _c, calculate and obtain by step [7]) and the imaging Output rusults in moment, u is input seismic wave field (be the conversion radio frequency channel collection data of input, obtained by step [1]).

In order to verify the effect of the inventive method, apply this inventive method actual 3D3C data is tested.

Data is from two bunch data of the three-dimensional three-component real data in new field, Sichuan.Fig. 6-1st, adopts the imaging road collection (CIP) of existing method direct rays anisotropy pre-stack time migration processing, and Fig. 6-2nd adopts the imaging road collection (CIP) of curved rays anisotropy pre-stack time migration of the present invention processing.Tu Zhong road collection is arranged (from left to right offset distance is descending) by far and near offset distance, and on CIP road, Fig. 6-1 collection, reflection line-ups upwarps gradually, and the inclined to one side phenomenon of obvious mistake (seeing the rectangle frame in Fig. 6-1) appears in road collection a long way.On Ze CIP road, Fig. 6-2 collection, there is obvious offset correction to even up effect for excessively inclined to one side a long way phenomenon.Two kinds of identical Migration velocity models of skew CIP road centralized procurement.The difference on migration imaging road collection that this causes just because of the whilst on tour error of calculation.The inventive method obtains obvious offset effect, has improved the image quality of offset gather far away.

Fig. 7 (a) is the imaging section that adopts the processing of existing method direct rays anisotropy pre-stack time migration, and Fig. 7 (b) is the imaging section that adopts patented method curved rays anisotropy pre-stack time migration of the present invention processing.On the migrated section figure of Fig. 7 (a), one group of reflection in 2.6 seconds can reliably be followed the trail of continuously hardly on section, and ripple group relation is indefinite.The reflected energy that in figure, square frame shows laterally suddenlys change, and low frequencyization is obvious.Adopt identical input data, treatment scheme, parameter and identical anisotropy Migration velocity model, on Fig. 7 (b) migrated section, within 2.6 seconds, one group of reflection wave groups relation is clear and definite, on section, can follow the trail of continuously explanation.The reflected energy that in figure, square frame shows is laterally even, and continuity is good, does not occur low frequency phenomenon, and the entire profile image quality is significantly improved.

Fig. 8 processes the collection comparison of CIP road without amplitude weight function and different weights function, and wherein, Fig. 8 (a) processes CIP road collection without weighting function, for the phenomenon that showed a long way weighting of large offseting distance.Fig. 8 (b) is that existing method simple distance amplitude weight function is processed CIP road collection, and to large offseting distance, weighting phenomenon appearred in processing a long way, close with the processing road collection without weighting function.Fig. 8 (c) is for adopting patented method weighting function formula of the present invention (12) to process CIP road collection, and amplitude variation with Offset meets AVO feature (seeing square frame in figure) more.

Fig. 9-1 and Fig. 9-2 are the Profile Correlation of three-dimensional PS ripple anisotropy pre-stack time migration processing, Fig. 9-1st, the pre-stack time migration section of the non-guarantor's width of existing method direct rays, Fig. 9-2nd, patented method curved rays of the present invention is protected the pre-stack time migration processing profiles of width.The pre-stack time migration processing profiles aggregate performance signal to noise ratio (S/N ratio) that curved rays is protected width is high, and deflection focusing energy is strong, migration is accurate, breakpoint is clear, substratum resolution is high.See figure interrupting layer and substratum explanation.

The present invention is that a 3D/2D transformed wave seismic data curved rays is protected width anisotropy prestack time migration method; the hypothesis of the method based on layer anisotropic media; in the Computation of seismic travel time formula of Kirchhoff migration; adopt 6 rank items and vertical axis of symmetry isotropy (VTI) parameter of offset distance; improve transformed wave whilst on tour computational accuracy, the migration processing of the multi-component seismic data that in adaptation, offset distance far away gathers.In amplitude of deflection weighting function calculates, according to the rate pattern estimation three-dimensional geometry invasin of time domain, the incident angle of accurate Calculation ray and emergence angle, estimation weighting factor carries out correction of amplitude to each imaging point, recovers seismic event real amplitude information.

Technique scheme is one embodiment of the present invention, for those skilled in the art, the invention discloses on the basis of application process and principle, be easy to make various types of improvement or distortion, and be not limited only to the described method of the above-mentioned embodiment of the present invention, therefore previously described mode is just preferred, and does not have restrictive meaning.

Claims (7)

1. transformed wave curved rays is protected a width anisotropy prestack time migration method, it is characterized in that: said method comprising the steps of:

(1) input three-dimensional transformed wave prestack road collection record and transformed wave anisotropy body of velocity, comprise that ripple vertical speed compares γ in length and breadth ₀, velocity equivalent compares γ _eff, transformed wave anisotropic parameters x _effwith transformed wave migration velocity V _cnthese four parameters;

(2) start a track data of input to carry out imaging processing;

(3) four parameters of applying step (1) input, are calculated to be the compressional wave migration velocity V at decent some place successively _pnwith shear wave migration velocity V _sn;

Calculate compressional wave anisotropic parameters η _effwith shear wave anisotropic parameters ξ _eff;

Calculate compressional wave whilst on tour t _p6with shear wave whilst on tour t _s6;

(4) the compressional wave curved rays whilst on tour t that applying step (3) is tried to achieve _p6, shear wave whilst on tour t _s6, compressional wave anisotropic parameters η _effwith shear wave anisotropic parameters ξ _effcalculate transformed wave anisotropy whilst on tour t _c;

(5) according to compressional wave migration velocity and compressional wave curved rays whilst on tour t _p6calculate ray incident angle α _scosine cos (α _s); According to shear wave migration velocity and shear wave curved rays whilst on tour t _s6calculate ray emergence angle α _gcosine cos (α _g);

(6) computational geometry diffusion term, L and sin (α);

(7) according to transformed wave amplitude weight function, calculated amplitude weighting coefficient W (ξ, R);

(8) according to transformed wave offset equation, application amplitude weight coefficient W (ξ, R) and transformed wave whilst on tour t _cand transformed wave input wave field derivative term, complete transformed wave curved rays and protect width migration processing;

(9) determine whether last imaging point, if so, proceed to step (10), if not, return to step (2);

(10) one imaging finish.

2. transformed wave curved rays according to claim 1 is protected width anisotropy prestack time migration method, it is characterized in that: in described step (3), be the compressional wave migration velocity V that is calculated to be respectively picture point place according to formula (4) _pnwith shear wave migration velocity V _sn:

(4)

In described step (3), be to be calculated to be respectively the compressional wave anisotropic parameters η of picture point place according to formula (5) _effwith shear wave anisotropic parameters ξ _eff:

(5)

In described step (3), be to calculate compressional wave whilst on tour t according to formula (10) and (11) _p6with shear wave whilst on tour t _s6:

Wherein, k is high-order small quantity, as constant coefficient item, and value 0.0～2.0;

T _p4for imaging point place compressional wave 4 rank item whilst on tours, t _s4for imaging point place shear wave 4 rank item whilst on tours;

Wherein, x _prepresent the distance from focal point to imaging point, x _srepresent the distance from geophone station to imaging point;

V in formula _i, h _kbe respectively interval velocity and zone thickness, n is hierarchy number.

3. transformed wave curved rays according to claim 2 is protected width anisotropy prestack time migration method, it is characterized in that: c _1p, c _2p, c _3pand c _4pfor the coefficient of compressional wave, adopt P-wave interval velocity to utilize formula (6) to calculate; c _1s, c _2s, c _3sand c _4sfor the coefficient of shear wave, adopt S-wave interval velocity to utilize formula (6) to calculate.

4. transformed wave curved rays according to claim 3 is protected width anisotropy prestack time migration method, it is characterized in that: described step (4) is to calculate transformed wave anisotropy whilst on tour t according to formula (12) _c:

Δ t in formula _pwith Δ t _scalculated by formula (2):

(2)。

5. transformed wave curved rays according to claim 4 is protected width anisotropy prestack time migration method, it is characterized in that: described step (5) utilizes formula (13) to calculate:

Tp in above formula, that Ts adopts is compressional wave whilst on tour t _p6with shear wave whilst on tour t _s6.

6. transformed wave curved rays according to claim 5 is protected width anisotropy prestack time migration method, it is characterized in that: described step (6) is achieved in that

First calculate offset midpoint position M (X _m, Y _m) coordinate:

(14)

X in formula _sand X _gbe respectively big gun inspection coordinate;

Then calculate reflection spot at ground projected position Q (X _q, Y _q) coordinate:

X _q＝CDP _x-X _m

Y _q＝CDP _y-Y _m (15)

CDP in formula _x, CDP _yfor the road header of present input data; L is the length of QM, is shown with end points coordinates table:

And the sine of the angle α of two straight line MQ and SG can calculate with the slope of two straight lines:

7. transformed wave curved rays according to claim 6 is protected width anisotropy prestack time migration method, it is characterized in that: described step (7) is to adopt formula (18) to calculate:

T in formula _rfor the time at imaging point place, T _pand T _s, the compressional wave whilst on tour t that adopts step (3) to calculate _p6with shear wave whilst on tour t _s6; V _cfor the migration velocity of imaging point, adopt described transformed wave migration velocity V _cn,, V _p0be compressional wave near-surface velocity, H is geophone offset.