CN104101901A - Converted-wave curved ray amplitude-reserved anisotropic pre-stack time offset time method - Google Patents
Converted-wave curved ray amplitude-reserved anisotropic pre-stack time offset time method Download PDFInfo
- Publication number
- CN104101901A CN104101901A CN201310116445.1A CN201310116445A CN104101901A CN 104101901 A CN104101901 A CN 104101901A CN 201310116445 A CN201310116445 A CN 201310116445A CN 104101901 A CN104101901 A CN 104101901A
- Authority
- CN
- China
- Prior art keywords
- wave
- whilst
- tour
- formula
- calculate
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Landscapes
- Geophysics And Detection Of Objects (AREA)
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
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 in formula
cfor C ripple whilst on tour, γ
0for compressional wave and shear wave vertical speed ratio, t
c0for C ripple t
0time, 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:
(2)
The form of weighting function is:
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
0, 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:
(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.
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 in formula
pwith Δ t
scalculated by formula (2):
(2)。
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.
Described step (6) is achieved in that
First calculate offset midpoint position M (Xm
,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:
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.
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 γ
0compare γ 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.
(4)
[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 γ
0compare γ 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.
(5)
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
1, c
2, c
3and c
4, 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,
In formula
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.
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.
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 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).
[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.):
(14)
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:
And the sine of the angle α of two straight line MQ and SG can calculate with the slope of two straight lines
[10] whole weighting function numerical procedure adopts the expression formula of formula (18),
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:
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
0, 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.
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201310116445.1A CN104101901A (en) | 2013-04-03 | 2013-04-03 | Converted-wave curved ray amplitude-reserved anisotropic pre-stack time offset time method |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201310116445.1A CN104101901A (en) | 2013-04-03 | 2013-04-03 | Converted-wave curved ray amplitude-reserved anisotropic pre-stack time offset time method |
Publications (1)
Publication Number | Publication Date |
---|---|
CN104101901A true CN104101901A (en) | 2014-10-15 |
Family
ID=51670194
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201310116445.1A Pending CN104101901A (en) | 2013-04-03 | 2013-04-03 | Converted-wave curved ray amplitude-reserved anisotropic pre-stack time offset time method |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN104101901A (en) |
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106154319A (en) * | 2015-04-22 | 2016-11-23 | 中国石油化工股份有限公司 | A kind of method for separating of imaging road collection |
CN109581494A (en) * | 2018-10-23 | 2019-04-05 | 中国石油天然气集团有限公司 | Prestack migration method and device |
CN112379431A (en) * | 2020-11-13 | 2021-02-19 | 中国地质科学院 | PS wave seismic data migration imaging method and system under complex surface condition |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20060072373A1 (en) * | 2004-10-04 | 2006-04-06 | Seongbok Lee | Amplitude preserving prestack migration method |
CN101285894A (en) * | 2008-05-30 | 2008-10-15 | 中国科学院地质与地球物理研究所 | Heaved earth surface collected seismic data direct prestack time migration method |
CN101315427A (en) * | 2007-05-29 | 2008-12-03 | 中国石油天然气集团公司 | Method and system for processing seismic exploration data of complex area |
-
2013
- 2013-04-03 CN CN201310116445.1A patent/CN104101901A/en active Pending
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20060072373A1 (en) * | 2004-10-04 | 2006-04-06 | Seongbok Lee | Amplitude preserving prestack migration method |
CN101315427A (en) * | 2007-05-29 | 2008-12-03 | 中国石油天然气集团公司 | Method and system for processing seismic exploration data of complex area |
CN101285894A (en) * | 2008-05-30 | 2008-10-15 | 中国科学院地质与地球物理研究所 | Heaved earth surface collected seismic data direct prestack time migration method |
Non-Patent Citations (5)
Title |
---|
CHUANWEN SUN ET AL.: "3D Kirchhoff PS-wave prestack time migration for V(z) and VTI media", 《EXPANDED ABSTRACTS OF 73RD ANNUAL INTERNAT SEG MTG》 * |
喻勤等: "基于CUDA的转换波Kirchhoff叠前时间偏移算法研究及实现", 《石油地球物理勘探》 * |
郝伟: "基于3_migs处理***的弯曲射线叠前时间偏移技术研究与应用", 《中国优秀硕士学位论文全文数据库 基础科学辑》 * |
陈枫: "基于GPU技术的叠前时间偏移及其在玛湖地区的应用", 《中国优秀硕士学位论文全文数据库 基础科学辑》 * |
黄中玉等: "三维多分量地震资料处理技术研究", 《石油物探》 * |
Cited By (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106154319A (en) * | 2015-04-22 | 2016-11-23 | 中国石油化工股份有限公司 | A kind of method for separating of imaging road collection |
CN106154319B (en) * | 2015-04-22 | 2018-05-11 | 中国石油化工股份有限公司 | A kind of method for separating for being imaged trace gather |
CN109581494A (en) * | 2018-10-23 | 2019-04-05 | 中国石油天然气集团有限公司 | Prestack migration method and device |
CN112379431A (en) * | 2020-11-13 | 2021-02-19 | 中国地质科学院 | PS wave seismic data migration imaging method and system under complex surface condition |
CN112379431B (en) * | 2020-11-13 | 2024-02-02 | 中国地质科学院 | PS wave seismic data migration imaging method and system under complex surface condition |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN101907727B (en) | Multi-component converted wave static correction method by using surface waves | |
CN104656142B (en) | One kind is using vertical seismic profiling (VSP) and the united seismic layer labeling method of well logging | |
CN108363101B (en) | A kind of inclined shaft crosshole seismic Gaussian beam pre-stack depth migration imaging method | |
CN101957455B (en) | Method of three-dimensional preserved-amplitude pre-stack time migration | |
CN102841375A (en) | Method for tomography velocity inversion based on angle domain common imaging gathers under complicated condition | |
CN102841379B (en) | Method for analyzing pre-stack time migration and speed based on common scatter point channel set | |
CN102213769A (en) | Method for determining anisotropic parameters by utilizing data of three-dimensional VSP (Vertical Seismic Profile) | |
CN105549081A (en) | Anisotropic medium common shot domain Gaussian beam migration imaging method | |
WO2012160409A1 (en) | Method of processing seismic data by providing surface offset common image gathers | |
CN104216014A (en) | Seismic signal frequency division processing method | |
CN102116869A (en) | High-precision prestack domain least square migration seismic imaging technology | |
Rosales et al. | Wave-equation angle-domain common-image gathers for converted waves | |
CN102053260B (en) | Method for acquiring azimuth velocity of primary wave and method for processing earthquake data | |
CN102636809A (en) | Method for generating spreading angle domain common image point gathers | |
CN105093301A (en) | Common imaging point reflection angle gather generation method and apparatus | |
CN102053262B (en) | Method for acquiring azimuth velocity of seismic converted wave and method for processing seismic data | |
CN102901984A (en) | Method for constructing true earth surface dip angle trace gathers of seismic data | |
CN104536041B (en) | Optimization method of seismological observation system parameters | |
CN103823241B (en) | A kind of method asking for reflecting Value of residual static correction in migration in offset domain | |
Zhu et al. | Recent applications of turning-ray tomography | |
CN104101901A (en) | Converted-wave curved ray amplitude-reserved anisotropic pre-stack time offset time method | |
CN102798888B (en) | Method for calculating velocity ratio of longitudinal wave to transverse wave by using non-zero wellhead distance data | |
Lafond et al. | Migration of wide‐aperture onshore‐offshore seismic data, central California: Seismic images of late stage subduction | |
Dhabaria et al. | Hamiltonian Monte Carlo based elastic full-waveform inversion of wide-angle seismic data | |
CN104076395A (en) | Mirror surface energy extraction and imaging method based on filtering combination |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
WD01 | Invention patent application deemed withdrawn after publication | ||
WD01 | Invention patent application deemed withdrawn after publication |
Application publication date: 20141015 |