CN106842300B - A kind of high efficiency multi-component seismic data true amplitude migration imaging method - Google Patents
A kind of high efficiency multi-component seismic data true amplitude migration imaging method Download PDFInfo
- Publication number
- CN106842300B CN106842300B CN201611197758.4A CN201611197758A CN106842300B CN 106842300 B CN106842300 B CN 106842300B CN 201611197758 A CN201611197758 A CN 201611197758A CN 106842300 B CN106842300 B CN 106842300B
- Authority
- CN
- China
- Prior art keywords
- wave field
- ray parameter
- time
- big gun
- parameter
- 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.)
- Expired - Fee Related
Links
- 238000013508 migration Methods 0.000 title claims abstract description 107
- 230000005012 migration Effects 0.000 title claims abstract description 107
- 238000003384 imaging method Methods 0.000 title claims abstract description 83
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 9
- 239000013598 vector Substances 0.000 claims description 31
- 238000004364 calculation method Methods 0.000 claims description 9
- 241000208340 Araliaceae Species 0.000 claims description 6
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 claims description 6
- 235000003140 Panax quinquefolius Nutrition 0.000 claims description 6
- 238000013461 design Methods 0.000 claims description 6
- 235000008434 ginseng Nutrition 0.000 claims description 6
- 238000005457 optimization Methods 0.000 claims description 6
- 238000005381 potential energy Methods 0.000 claims description 6
- 238000005070 sampling Methods 0.000 claims description 6
- 238000007689 inspection Methods 0.000 claims description 5
- 238000004321 preservation Methods 0.000 claims description 4
- 238000006073 displacement reaction Methods 0.000 claims description 3
- 238000013213 extrapolation Methods 0.000 claims description 3
- 230000002035 prolonged effect Effects 0.000 claims 1
- 238000000034 method Methods 0.000 abstract description 44
- 238000012545 processing Methods 0.000 abstract description 5
- 238000005516 engineering process Methods 0.000 abstract description 2
- 230000011218 segmentation Effects 0.000 abstract description 2
- 238000007796 conventional method Methods 0.000 description 10
- 230000000694 effects Effects 0.000 description 6
- 238000004880 explosion Methods 0.000 description 3
- 150000003839 salts Chemical class 0.000 description 3
- 230000001066 destructive effect Effects 0.000 description 2
- 238000012795 verification Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- BULVZWIRKLYCBC-UHFFFAOYSA-N phorate Chemical compound CCOP(=S)(OCC)SCSCC BULVZWIRKLYCBC-UHFFFAOYSA-N 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 230000005477 standard model Effects 0.000 description 1
- 239000011800 void material Substances 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
A kind of high efficiency multi-component seismic data true amplitude migration imaging method, applied to multi-component earthquake data processing technology field, this method is improved on the basis of elastic wave reverse-time migration method, it is characterized in being directly input with multi-component seismic data, pass through the measures such as encryption algorithm, segmentation imaging, again under the frame of inverting, the elastic wave reverse-time migration based on inverting is realized.Coding, time slice and the thought of inverting are introduced into elastic wave reverse-time migration by the present invention, and compared with common elastic wave reverse-time migration, the present invention can obtain the pre-stack depth migration section of high-precision and hi-fi of amplitude, while the efficiency of algorithm is very high.The present invention completely keeps vectorial property, amplitude and the phase characteristic of longitudinal and shear wave, migrated section to have fidelity；It effectively eliminates and deviates illusion caused by crosstalk between longitudinal and shear wave, improve the precision of imaging；The amount of storage for effectively reducing elastic wave reverse-time migration, improves computational efficiency.
Description
Technical field
The invention belongs to field of geophysical exploration, are related to the processing of multi-component seismic data prestack migration image, especially relate to
And a kind of high efficiency, low storage multi-component seismic data true amplitude migration imaging method.
Background technology
Based on longitudinal wave exploration method, this method thinks to only exist in underground medium vertical the overwhelming majority in traditional oil-gas exploration
Wave, however, practically descending the wave field of Propagation not only to have longitudinal wave, but also the also wave modes such as shear wave and converted wave, i.e.,
Seismic wave field is elastic wave field.Multi-component seismic data contains longitudinal wave naturally, and a variety of wave modes such as shear wave and converted wave use
Traditional simple component seismic data processing technique certainly exists larger error to handle multi-component seismic data, even results in processing
As a result mistake.Pre-stack depth migration is a vital link in seismic processing chain.Current existing multi -components
Seismic data pre-stack depth migration imaging method can be divided into three classes：Ray class offset method, one-way wave class offset method and bullet
Property wave reverse-time migration method.Wherein, elastic wave reverse-time migration is precision highest, the most steady one kind of algorithm.Nevertheless, elastic
Wave reverse-time migration still remains many disadvantages, such as：Longitudinal and shear wave crosstalk causes serious offset illusion, the reversion of shear wave polarity to cause
Migrated section there are destructive interference, migrated section amplitude infidelity, migrated section resolution ratio is low, migration algorithm computational efficiency is low
Deng if in addition, using wave separation technique in length and breadth, the amplitude of wave field and phase information being caused to change, change wave field
Kinetic characteristics.The prestack depth that these described problems make elastic wave reverse-time migration more difficult applied to practical multi-component seismic data
In degree offset.For this purpose, a set of new, with the characteristics of " high efficiency, protects width at high-precision " high efficiency multi-component seismic money must be established
Expect real amplitude pre-stack depth migration imaging method.
Invention content
It is an object of the invention in view of the foregoing defects the prior art has, provide it is a kind of it is based on equations for elastic waves,
Multi-component seismic data pre-stack depth migration imaging method with the characteristics of " high-precision, low storage, protects width at high efficiency ".
In order to achieve the goal above, the present invention provides a kind of high efficiency multi-component seismic data true amplitude migration imaging side
Method includes the following steps：
(a)：The more big gun multi component seismic records D (x obtained by field inspectionr, t；xs)=(Dx, Dy, Dz), give background
Model, parameter preset determine the multi component seismic records for carrying out migration imaging, observation system parameter and source function S (x, t；
xs), wherein x representation space position vectors, xsIndicate hypocentral location vector, xrIndicate that geophone station position vector, t indicate that the time becomes
Amount；
(b)：Observation system parameter based on reading, design plane wave coding function C (xs, p), specially：
C(xs, p) and=fδ(t-p·xs), (1)
Wherein, p is ray parameter variable, fδIndicate impulse function；Plane wave coding function C (x hereins, p) and when being equivalent to
Move function, time shift amount pxsWith focal point position xsDifferent and change, ray parameter p is specially：
Wherein, α is earth's surface velocity of longitudinal wave, and θ is incidence angle；By different incidence angle θs, different ray ginsengs can be obtained
Number p；
By the coding function, and using each ray parameter p to more big gun multi component seismic records D (x of observationr,
t；xs) encoded, obtain the coding multi -components big gun collection D about each ray parameter pC(xr,t；P), specially：
DC(xr,t；P)=∫ D (xr,t；xs)·C(xs,p)dxs； (3)
Equally, using identical coding function, and using each ray parameter p to observation system described in step (a)
In source function S (x, t；xs) encoded, obtain the coded source function S of each ray parameter pC(x,t；P), specially：
SC(x,t；P)=∫ S (x, t；xs)·C(xs,p)dxs； (4)
(c)：Utilize background model and the dominant record time t of more big gun multi component seismic records of observationmax, design record
The piecewise interval of time t, specially：It will record time t point be three sections, be indicated with following formula：
Wherein, tR1, tR2And tR3Three interval ranges of record time, t are indicated respectively1(x) and t2(x) specific calculating formula
For：
Wherein, αmaxFor the maximum speed of background model, hzRepresent depth, hxRepresent horizontal distance；
(d)：For each ray parameter p, it is performed both by following operation：Using each ray parameter p, pass through step
(b) coding mode in obtains coded source function corresponding with the ray parameter, the bullet based on conventional isotropic medium
Property wave equation obtains the ray parameter pair using background model given herein above to the coded source into traveling-wave field forward direction continuation
The source wavefield at each moment answered, in each time step of the source wavefield continuation, when using record in step (c)
Between piecewise interval as time constraint condition, when each using conventional elastic wave field energy balane equation calculation source wavefield
It carves, the wavefield energy value of each position, is judged by energy values size, obtained in each time interval, each position wave field
At the time of corresponding to Energy maximum value, using this moment as the imaging time T (x) in the time interval, corresponding to the position,
At the moment, the position, using image-forming condition formula, preserve with the relevant item of source wavefield in image-forming condition, while preserving should
Moment, the position wave field maximum energy value；
And then obtain in the corresponding imaging time of each ray parameter, image-forming condition with the relevant item of source wavefield and wave
Field maximum energy value；Finally the wave field maximum energy value in each corresponding period of each ray parameter is folded
Add, obtains illuminance compensation factor；
(e)：For each ray parameter p, it is performed both by following operation：It is same based on described to each ray parameter p
Isotropic medium equations for elastic waves, to step (b) the coding multi -components big gun collection D that ray parameter p is correspondingC(xr,
t；P) into the continuation of traveling-wave field inverse time, the geophone station wave field at each corresponding moment of the ray parameter is obtained, in the geophone station wave
In each time step of continuation, using the imaging time T (x) of each period, each position that are preserved in step (d),
Judge whether current time meets imaging time, if conditions are not met, then without imaging, if it is satisfied, then according to image-forming condition,
It is imaged, is somebody's turn to do with the relevant item of source wavefield and geophone station wave field using in the image-forming condition preserved in step (d)
The corresponding migrated section of ray parameter；
And then obtain the corresponding migrated section of each ray parameter；By all corresponding migrated sections of ray parameter into
Row superposition, while energy compensating is carried out to the migrated section of superposition using illuminance compensation factor in step (d), it obtains initial inclined
Move section namely the migrated section of the 0th iteration
(f)：For each ray parameter p, following operation is executed：To each ray parameter p, by step (b)
Coding mode obtain corresponding with ray parameter coded source function, the equally elasticity based on the isotropic medium
Wave equation obtains ray parameter correspondence to the coded source using background model given herein above into traveling-wave field forward direction continuation
Background wave field；Utilize the initial offset section I described in step (e)0, the ray parameter background wave field and step
Suddenly the background model given in (a), builds the corresponding virtual focus of the ray parameter；Based on the virtual focus it is each to
Same sex dielectric resilient wave equation obtains the ray parameter to the shot point using given background model into traveling-wave field forward direction continuation
Corresponding inverse migration wave field carries out record sampling to the corresponding inverse migration wave field of the ray parameter, obtains ray parameter correspondence
Coding multi -components inverse migration big gun collection；
And then obtain the corresponding coding multi -components inverse migration big gun collection of each ray parameter, the coding multi -components big gun of all big guns
Collection constitutes initial code multi -components inverse migration big gun collection namely the coding multi -components inverse migration big gun collection of the 0th iterationAt this point, it is i=0 that current iterations, which are arranged,；
(g)：Current iterations i=i+1 is set, based on the migrated section obtained by (i-1)-th iterationWith the coding multi -components inverse migration big gun collection obtained by (i-1)-th iterationCarry out inverting iteration；
It is for each ray parameter p, the ray parameter obtained by (i-1)-th iteration is corresponding in ith iteration
Encode multi -components inverse migration big gun collection di-1The coding multi -components big gun collection D of observation corresponding with the ray parameterCIt is poor to make, and calculates this and penetrates
The corresponding coding multi -components big gun collection residual error of line parameter；And then obtain the corresponding coding multi -components big gun collection residual error of all ray parameters
It is described by the corresponding coding multi -components inverse migration big gun collection d of the ray parameter obtained by (i-1)-th iterationi-1With this
The coding multi -components big gun collection D of the corresponding observation of ray parameterCIt is poor to make, and it is residual to calculate the corresponding coding multi -components big gun collection of the ray parameter
Difference, specially：
Δdi-1=di-1-DC； (7)
(h)：For each ray parameter p, following operation is executed：To each ray parameter p, equally based on described
The equations for elastic waves of isotropic medium, coding multi -components big gun collection residual error corresponding to the ray parameterInto the continuation of traveling-wave field inverse time, each corresponding moment of the ray parameter is obtained
Geophone station wave field, in each time step of the geophone station wave field extrapolation, using each period preserved in step (d),
The imaging time T (x) of each position, judges whether current time meets imaging time, if conditions are not met, then without being imaged,
If it is satisfied, then according to image-forming condition, using in the image-forming condition preserved in step (d) with the relevant item of source wavefield and inspection
Wave point wave field is imaged, and the corresponding gradient profile of the ray parameter is obtained；
And then obtain the corresponding gradient profile of each ray parameter；By all corresponding gradient profiles of ray parameter into
Row superposition, while energy compensating is carried out to the gradient profile of superposition using illuminance compensation factor in step (d), it obtains this and changes
The gradient profile in generation namely the gradient profile of ith iteration
(i)：Using optimal inversion algorithm, the descent direction section of ith iteration is built
(j)：For each ray parameter p, obtained by the coding mode in step (b) corresponding with the ray parameter
Coded source function, the equally equations for elastic waves based on the isotropic medium, utilize background model pair given herein above
The coded source obtains the corresponding background wave field of the ray parameter into traveling-wave field forward direction continuation；Obtained in step (i)
The descent direction section r of ith iterationi, the ray parameter background wave field and step (a) in give background mould
Type builds the corresponding virtual focus of the ray parameter；Equally based on the isotropic medium elastic wave with the virtual focus
It is reverse-biased to obtain the corresponding disturbance of the ray parameter into traveling-wave field forward direction continuation to the shot point using given background model for equation
Wave field is moved, disturbance inverse migration wave field corresponding to the ray parameter carries out record sampling, obtains the corresponding coding of the ray parameter
Multi -components disturb inverse migration big gun collection；
And then obtain the corresponding coding multi -components disturbance inverse migration big gun collection of each ray parameter, more points of the coding of all big guns
Amount disturbance big gun collection constitutes the coding multi -components disturbance inverse migration big gun collection of ith iteration
Stepsize formula is recycled to calculate the optimization step-length α of current iterationi；
The stepsize formula is
In the equation (8), summation variable k includes tri- directions x, y and z of cartesian coordinate；
(k)：The optimization step-length α obtained using step (j)iAnd the descent direction section r that step (i) obtainsi, update ith
The migrated section I of iterationi=Ii-1+αiri；Optimize step-length α using statingiAnd the coding multi -components disturbance inverse migration that step (j) obtains
Big gun collection δ di, update the coding multi -components inverse migration big gun collection d of ith iterationi=di-1+αiδdi；
The migrated section I of the update ith iterationi=Ii-1+αiriSpecially：
The multi -components inverse migration of the update ith iteration records di=di-1+αiδdiSpecially：
(l)：After ith iteration is complete, the cost functional value f of ith iteration is calculatedi, judge whether current iteration meets receipts
Standard is held back, it is final migrated section that newest migrated section is exported if meetingOtherwise step (g) is repeated
To (k), until obtaining final migrated section；
The cost functional value f for calculating ith iterationiSpecially：
In the equation (11), summation variable k includes tri- directions x, y and z of cartesian coordinate；
The convergence is specially：
In the equation (12), theshold indicates the threshold value standard of iteration stopping, usually chooses 0.001.
The present invention implement use technical solution further include：In step (d), step (e), step (f), step (h) and step
Suddenly the equations for elastic waves of the isotropic medium described in (j) is：
In the equation (13), V=(Vx,Vy,Vz)TIndicate vector velocity wave field；σ=(σxx,σyy,σzz,σxy,σyz,σxz)T
Indicate vector stress wave field；λ and μ indicates medium Lame constants；ρ indicates Media density；X, y and z indicate cartesian coordinate respectively
Tri- directions x, y and z.
The present invention implement use technical solution further include：In step (d), the imaging described in step (e) and step (h)
Condition is：
Wherein, α is medium velocity of longitudinal wave, and β is medium shear wave velocity, and the p-and s-wave velocity of medium and the medium Lame are normal
Meet relationship between number, Media density：
With
In the equation (14), in step (e),In step (h),
In the equation (11),Indicate the vector velocity wave field in source wavefield,Indicate the vector stress wave field in geophone station wave field.
The present invention implement use technical solution further include：Described in step (d) in each of the source wavefield continuation
In a time step, using the record time slice section in step (c) as time constraint condition, elastic wave field energy meter is utilized
Calculate equation calculation each moment, each position wavefield energy value, specially：Elastic wave field energy balane equation is：
E=K+W, (15)
Wherein, E is the energy of elastic wave field, and K is the kinetic energy of elastic wave field, and W is the potential energy of elastic wave field；Elastic wave field
Kinetic energy and potential energy can be obtained by its velocity field and stress field calculation, specially：
Wherein, U is displacement field, which can be obtained by doing time integral to the vector velocity wave field V；
Judged by energy values size described in step (d), is obtained in each time interval, each position wave field energy
At the time of measuring corresponding to maximum value, using this moment as the imaging time T (x) in the time interval, corresponding to the position, specifically
For：With first time interval t in step (c)R1For, in each position, if current time, t was in tR1Within, then
The wavefield energy E at current timet(x) if more than in the time interval, the maximum value E of the wavefield energy at all moment beforemax
(x), then, Emax(x)=Et(x), the imaging time T (x) of the time interval=t；Otherwise, in the time interval, the position
The maximum value and imaging time of wavefield energy remain unchanged；For other two time interval, then similarly；
Described in step (d) at the moment, the position, using image-forming condition formula, preserve in image-forming condition with shake
The relevant item of source wave field, at the same preserve the moment, the position wave field maximum energy value, specially：In image-forming condition with focus
The relevant item of wave field includes： WithTotally nine, meter
It calculates the moment, this nine value and preserve the value of calculating at the position；Meanwhile preserve the moment, the position wave field
Maximum energy value Emax(x)。
The present invention implement use technical solution further include：Utilization step (d) described in step (e) and step (h)
Each period of middle preservation, the imaging time T (x) of each position, judge whether current time meets imaging time, if not
Meet, then without imaging, if it is satisfied, then according to image-forming condition, using in the image-forming condition preserved in step (d) with focus
The relevant item of wave field and geophone station wave field are imaged, and obtain the corresponding migrated section of the ray parameter, specially：With current
Moment t meets first time interval t in step (c)R1For, in each position, if current time t satisfaction：|T(x)-t
| < 5, then according to the image-forming condition, to the phase of the geophone station wave field obtained by the source wavefield continuous item of preservation and current time
It closes item and applies the image-forming condition, obtain m1And m2, and the imaging results at the moment are added in the imaging section of the position；
Otherwise, the moment, the position are without imaging；For other two time interval, then similarly.
The present invention implement use technical solution further include：Based on the void described in step (f) and step (j)
The isotropic medium equations for elastic waves of quasi- focus, using given background model to the shot point into traveling-wave field forward direction continuation, this
The isotropic medium equations for elastic waves with the virtual focus at place is expressed as：
In the equation (17), δ V=(δ Vx,δVy,δVz)TIndicate inverse migration vector velocity wave field；δ σ=(δ σxx,δσyy,δ
σzz,δσxy,δσyz,δσzx)TIndicate inverse migration vector stress wave field；
In the equation (17), virtual source vector F=(Fxx,Fyy,Fzz,Fxy,Fyz,Fxz)TIt is expressed as：
In the equation (18), intermediate variable m '1With m '2With the m1And m2Meet following relationship：
With
In step (f),In step (j),
The invention adopts the above technical scheme, which has the following advantages：1) thought of inverting is introduced bullet by the present invention
In property wave reverse-time migration, compared with common elastic wave reverse-time migration, the present invention can obtain high-precision, high-resolution, high noise
Than the migrated section of, hi-fi of amplitude；2) image-forming condition used in the present invention corrects shear wave polarity naturally, and school is inverted without polarity
Migrated section destructive interference caused by positive operation can effectively overcome shear wave polarity to invert；3) imaging mode that the present invention uses,
Operation is detached without additional longitudinal and shear wave, and completely maintains vectorial property, amplitude and the phase characteristic of longitudinal and shear wave, is had
It eliminates to effect between longitudinal and shear wave and deviates illusion caused by crosstalk, greatly improve the precision of imaging；4) present invention gained at
Picture section physical significance is very clear, is convenient for later stage GEOLOGICAL INTERPRETATION；5) it present invention can be widely used in petroleum exploration domain, it is special
It is not more obvious for the imaging effect in complicated structure deep.
Description of the drawings
Fig. 1 is high efficiency multi-component seismic data true amplitude migration imaging method flow diagram provided by the invention；
Fig. 2 is the hollow illustraton of model of two dimension provided by the invention；Wherein, Fig. 2 (a) is velocity of longitudinal wave model；Fig. 2 (b) is horizontal
Velocity model；Fig. 2 (c) is background velocity of longitudinal wave model；Fig. 2 (d) is background shear wave velocity model；
Fig. 3 is the coding big gun collection after excision direct wave provided by the invention, corresponding ray parameter p=0 μ s/m：Wherein,
Fig. 3 (a) is the observational record horizontal component of coding；Fig. 3 (b) is the observational record vertical component of coding；Fig. 3 (c) is to utilize this
The initial inverse migration recording level component of coding obtained by inventive method；Fig. 3 (d) is to utilize the coding obtained by the method for the present invention
Initial inverse migration records vertical component；Fig. 3 (e) is the final inverse migration recording level point using the coding obtained by the method for the present invention
Amount；Fig. 3 (f) is the final inverse migration record vertical component using the coding obtained by the method for the present invention；
Fig. 4 is more big guns superposition migrated section of the hollow model of two dimension shown in Fig. 2：Wherein, Fig. 4 (a) is to utilize conventional method
The I of gainedαSection, Fig. 4 (b) are to utilize the I obtained by conventional methodβSection, Fig. 4 (c) are to utilize the I obtained by the method for the present inventionαIt cuts open
Face, Fig. 4 (d) are to utilize the I obtained by the method for the present inventionβSection；
Fig. 5 is Marmousi-ii models provided by the invention：Wherein, Fig. 5 (a) is velocity of longitudinal wave model；Fig. 5 (b) is horizontal
Velocity model；Fig. 5 (c) is background velocity of longitudinal wave model；Fig. 5 (d) is background shear wave velocity model；
Fig. 6 is more big guns superposition migrated section of Marmousi-ii models shown in Fig. 5：Wherein, Fig. 6 (a) is to utilize tradition side
I obtained by methodαSection；Fig. 6 (b) is to utilize the I obtained by conventional methodβSection；Fig. 6 (c) is to utilize the I obtained by the method for the present inventionα
Section；Fig. 6 (d) is to utilize the I obtained by the method for the present inventionβSection；
Fig. 7 is SEG/EAGE Salt models provided by the invention：Wherein, Fig. 7 (a) is velocity of longitudinal wave model；Fig. 7 (b) is
Shear wave velocity model；Fig. 7 (c) is background velocity of longitudinal wave model；Fig. 7 (d) is background shear wave velocity model；
Fig. 8 is more big guns superposition migrated section of the models of SEG/EAGE Salt shown in Fig. 7：Wherein, Fig. 8 (a) is to utilize tradition
I obtained by methodαSection；Fig. 8 (b) is to utilize the I obtained by conventional methodβSection；Fig. 8 (c) is using obtained by the method for the present invention
IαSection；Fig. 8 (d) is to utilize the I obtained by the method for the present inventionβSection；
Specific implementation mode
In order to the purpose of the present invention, technical solution and advantage are more clearly understood, and once combine accompanying drawings and embodiments, right
The present invention is further elaborated.It should be appreciated that specific embodiment described herein is only used for explaining the present invention, not
For limiting the present invention.
Referring to Fig. 1, being that high efficiency multi-component seismic data true amplitude migration imaging method flow provided by the invention is shown
It is intended to.The high efficiency multi-component seismic data true amplitude migration imaging method of the embodiment of the present invention includes the following steps：
Step S100：The more big gun multi component seismic records D (x obtained by field inspectionr,t；xs)=(Dx,Dy,Dz), it gives
Determine background model, parameter preset, determines multi component seismic records, observation system parameter and the source function S for carrying out migration imaging
(x,t；xs), wherein x representation space position vectors, xsIndicate hypocentral location vector, xrIndicate that geophone station position vector, t indicate
Time variable；
Step S200：Observation system parameter based on reading, design plane wave coding function C (xs, p), specially：
C(xs, p) and=fδ(t-p·xs), (1)
Wherein, p is ray parameter variable, fδIndicate impulse function；Plane wave coding function C (x hereins, p) and when being equivalent to
Move function, time shift amount pxsWith focal point position xsDifferent and change, ray parameter p is specially：
Wherein, α is earth's surface velocity of longitudinal wave, and θ is incidence angle；By different incidence angle θs, different ray ginsengs can be obtained
Number p；
By the coding function, and using each ray parameter p to more big gun multi component seismic records D (x of observationr,
t；xs) encoded, obtain the coding multi -components big gun collection D about each ray parameter pC(xr,t；P), specially：
DC(xr,t；P)=∫ D (xr,t；xs)·C(xs,p)dxs； (3)
Equally, using identical coding function, and using each ray parameter p to observation system described in step (a)
In source function S (x, t；xs) encoded, obtain the coded source function S of each ray parameter pC(x,t；P), specially：
SC(x,t；P)=∫ S (x, t；xs)·C(xs,p)dxs； (4)
Step S300：Utilize background model and the dominant record time t of more big gun multi component seismic records of observationmax, design
The piecewise interval of time t is recorded, specially：It will record time t point be three sections, be indicated with following formula：
Wherein, tR1, tR2And tR3Three interval ranges of record time, t are indicated respectively1(x) and t2(x) specific calculating formula
For：
Wherein, αmaxFor the maximum speed of background model, hzRepresent depth, hxRepresent horizontal distance；
Step S400：For each ray parameter p, it is performed both by following operation：Using each ray parameter p, pass through
Coding mode in step S200 obtains coded source function corresponding with the ray parameter, based on conventional isotropic medium
Equations for elastic waves obtain ray ginseng using background model given herein above to the coded source into traveling-wave field forward direction continuation
The source wavefield at number corresponding each moment, in each time step of the source wavefield continuation, using in step S300
Time slice section is recorded as time constraint condition, the elastic wave field energy balane equation calculation source wavefield using routine is every
A moment, each position wavefield energy value, judged by energy values size, obtained in each time interval, each position
At the time of corresponding to wavefield energy maximum value, using this moment as the imaging time T in the time interval, corresponding to the position
(x), at the moment, the position, using image-forming condition formula, preserve in image-forming condition with the relevant item of source wavefield, simultaneously
Preserve the moment, the position wave field maximum energy value；And then obtain the corresponding imaging time of each ray parameter, at slice
In part with the relevant item of source wavefield and wave field maximum energy value；Finally by each corresponding period of each ray parameter
Interior wave field maximum energy value is overlapped, and obtains illuminance compensation factor；
The equations for elastic waves of conventional isotropic medium described in step S400 is：
In the equation (9), V=(Vx,Vy,Vz)TIndicate vector velocity wave field；σ=(σxx,σyy,σzz,σxy,σyz,σxz)T
Indicate vector stress wave field；λ and μ indicates medium Lame constants；ρ indicates Media density；X, y and z indicate cartesian coordinate respectively
Tri- directions x, y and z；
Image-forming condition described in step S400 is：
Wherein, α is medium velocity of longitudinal wave, and β is medium shear wave velocity, and the p-and s-wave velocity of medium and the medium Lame are normal
Meet relationship between number, Media density：
With
Described in step S400 in each time step of the source wavefield continuation, utilize the note in step S300
Time slice section is recorded as time constraint condition, the elastic wave field energy balane equation calculation source wavefield using routine is each
Moment, each position wavefield energy value, specially：Elastic wave field energy balane equation is：
E=K+W, (15)
Wherein, E is the energy of elastic wave field, and K is the kinetic energy of elastic wave field, and W is the potential energy of elastic wave field；Elastic wave field
Kinetic energy and potential energy can be obtained by its velocity field and stress field calculation, specially：
Wherein, U is displacement field, which can be obtained by doing time integral to the vector velocity wave field V；
Described is judged by energy values size, is obtained in each time interval, each position wavefield energy maximum value
At the time of corresponding, using this moment as the imaging time T (x) in the time interval, corresponding to the position, specially：With step
First time interval t in rapid S300R1For, in each position, if current time, t was in tR1Within, then when current
The wavefield energy E at quartert(x) if more than in the time interval, the maximum value E of the wavefield energy at all moment beforemax(x), that
, Emax(x)=Et(x), the imaging time T (x) of the time interval=t；Otherwise, in the time interval, the wave field energy of the position
The maximum value and imaging time of amount remain unchanged；For other two time interval, then similarly；
It is described at the moment, the position, using image-forming condition formula, preserve related to source wavefield in image-forming condition
Item, while preserve the moment, the position wave field maximum energy value, specially：It is relevant with source wavefield in image-forming condition
Include： WithTotally nine, calculate the moment,
This nine value and the value of calculating is preserved at the position；Meanwhile preserve the moment, the position wave field maximum energy value
Emax(x)；
Step S500：For each ray parameter p, it is performed both by following operation：To each ray parameter p, same base
In the equations for elastic waves of the isotropic medium, the coding multi -components big gun collection D corresponding to ray parameter pC(xr,t；
P) into the continuation of traveling-wave field inverse time, the geophone station wave field at each corresponding moment of the ray parameter is obtained, in the geophone station wave field
In each time step of continuation, using the imaging time T (x) of each period, each position that are preserved in step S400, sentence
Whether disconnected current time meets imaging time, if conditions are not met, then without imaging, if it is satisfied, then according to image-forming condition, profit
It is imaged with the relevant item of source wavefield and geophone station wave field in the image-forming condition preserved in step S400, obtains this and penetrate
The corresponding migrated section of line parameter；And then obtain the corresponding migrated section of each ray parameter；By all ray parameters pair
The migrated section answered is overlapped, while carrying out energy to the migrated section of superposition using illuminance compensation factor in step S400
Compensation obtains initial offset section namely the migrated section of the 0th iteration
The imaging time T using each period, each position that are preserved in step S400 described in step S500
(x), judge whether current time meets imaging time, if conditions are not met, then without imaging, if it is satisfied, then according to imaging
Condition is imaged using in the image-forming condition preserved in step S400 with the relevant item of source wavefield and geophone station wave field,
The corresponding migrated section of the ray parameter is obtained, specially：Meet first time interval t in step S300 with current time tR1
For, in each position, if current time t satisfaction：| T (x)-t | < 5, then according to the image-forming condition, the shake to preservation
The continuous item of geophone station wave field obtained by source wave field continuous item and current time applies the image-forming condition, obtains m1And m2, and will
The imaging results at the moment are added in the imaging section of the position；Otherwise, the moment, the position are without imaging；For another
Outer two time intervals, then similarly；
Image-forming condition is applied described in step S500, is specially based on equation (14), and migrated section is specially：
Step S600：For each ray parameter p, following operation is executed：To each ray parameter p, pass through step
Coding mode in S200 obtains coded source function corresponding with the ray parameter, equally based on the isotropic medium
Equations for elastic waves obtain ray ginseng using background model given herein above to the coded source into traveling-wave field forward direction continuation
The corresponding background wave field of number；Utilize the initial offset section I described in step S5000, the ray parameter background wave
And step S100 in give background model, build the corresponding virtual focus of the ray parameter；Based on described virtual
The isotropic medium equations for elastic waves of focus obtains the shot point into traveling-wave field forward direction continuation using given background model
The corresponding inverse migration wave field of the ray parameter carries out record sampling to the corresponding inverse migration wave field of the ray parameter, obtains this and penetrate
The corresponding coding multi -components inverse migration big gun collection of line parameter；And then obtain the corresponding coding multi -components inverse migration of each ray parameter
The coding multi -components big gun collection of big gun collection, all big guns constitutes initial code multi -components inverse migration big gun collection namely the coding of the 0th iteration
Multi -components inverse migration big gun collectionAt this point, it is i=0 that current iterations, which are arranged,；
The isotropic medium equations for elastic waves with virtual focus described in step S600 is expressed as：
In the equation (17), δ V=(δ Vx,δVy,δVz)TIndicate inverse migration vector velocity wave field；δ σ=(δ σxx,δσyy,δ
σzz,δσxy,δσyz,δσzx)TIndicate inverse migration vector stress wave field；In the equation (17), virtual source vector F=(Fxx,
Fyy,Fzz,Fxy,Fyz,Fxz)TIt is expressed as：
In the equation (18), intermediate variable m '1With m '2With the m1And m2Meet following relationship：
With
In step S600, in the equation (18),
Step S700：Current iterations i=i+1 is set, based on the migrated section obtained by (i-1)-th iterationWith the coding multi -components inverse migration big gun collection obtained by (i-1)-th iterationCarry out inverting iteration；
It is for each ray parameter, the ray parameter obtained by (i-1)-th iteration is corresponding in ith iteration
Encode multi -components inverse migration big gun collection di-1The coding multi -components big gun collection D of observation corresponding with the ray parameterCIt is poor to make, and calculates this and penetrates
The corresponding coding multi -components big gun collection residual error of line parameter；And then obtain the corresponding coding multi -components big gun collection residual error of all ray parameters
It is described by the corresponding coding multi -components inverse migration big gun collection d of the ray parameter obtained by (i-1)-th iterationi-1With this
The coding multi -components big gun collection D of the corresponding observation of ray parameterCIt is poor to make, and it is residual to calculate the corresponding coding multi -components big gun collection of the ray parameter
Difference, specially：
Δdi-1=di-1-DC； (7)
Step S800：For each ray parameter p, following operation is executed：To each ray parameter p, equally it is based on
The equations for elastic waves of the isotropic medium, coding multi -components big gun collection residual error corresponding to the ray parameterInto the continuation of traveling-wave field inverse time, each corresponding moment of the ray parameter is obtained
Geophone station wave field, in each time step of the geophone station wave field extrapolation, using each period preserved in step S400,
The imaging time T (x) of each position, judges whether current time meets imaging time, if conditions are not met, then without being imaged,
If it is satisfied, then according to image-forming condition, using in the image-forming condition preserved in step S400 with the relevant item of source wavefield and
Geophone station wave field is imaged, and the corresponding gradient profile of the ray parameter is obtained；And then it is corresponding to obtain each ray parameter
Gradient profile；All corresponding gradient profiles of ray parameter are overlapped, while being compensated using illuminance in step S400
The gradient profile of factor pair superposition carries out energy compensating, and the gradient of the gradient profile namely ith iteration that obtain current iteration is cutd open
Face
It is described to apply image-forming condition in step S800, specially it is based on the equation (14), wherein
Step S900：Using optimal inversion algorithm, the descent direction section of ith iteration is built
Step S1000：For each ray parameter p, is obtained by the coding mode in step S200 and joined with the ray
The corresponding coded source function of number, the equally equations for elastic waves based on the isotropic medium, utilize the back of the body given herein above
Scape model, into traveling-wave field forward direction continuation, obtains the corresponding background wave field of the ray parameter to the coded source；Utilize step S900
Obtained in ith iteration descent direction section ri, the ray parameter background wave field and step S100 in
Given background model builds the corresponding virtual focus of the ray parameter；Equally based on each to same of the virtual focus
Property dielectric resilient wave equation obtains the ray parameter pair using given background model to the shot point into traveling-wave field forward direction continuation
The disturbance inverse migration wave field answered, disturbance inverse migration wave field corresponding to the ray parameter carry out record sampling, obtain ray ginseng
The corresponding coding multi -components of number disturb inverse migration big gun collection；And then it is anti-to obtain the corresponding coding multi -components disturbance of each ray parameter
Big gun collection is deviated, the coding multi -components disturbance big gun collection of all big guns constitutes the coding multi -components disturbance inverse migration big gun collection of ith iterationStepsize formula is recycled to calculate the optimization step-length α of current iterationi；
The descent direction section r using the ith iteration obtained in step S900 described in step S1000i, institute
It is corresponding virtual to build the ray parameter for the background model given in the background wave field and step S100 of the ray parameter stated
Focus is specially based on the equation (17), wherein
Stepsize formula described in step S1000 is
In the equation (8), summation variable k includes tri- directions x, y and z of cartesian coordinate；
Step S1100：The optimization step-length α obtained using step S1000iAnd the descent direction section r that step S900 is obtainedi,
Update the migrated section I of ith iterationi=Ii-1+αiri；Optimize step-length α using statingiAnd more points of the coding that step S1000 is obtained
Amount disturbance inverse migration big gun collection δ di, update the coding multi -components inverse migration big gun collection d of ith iterationi=di-1+αiδdi；
The migrated section I of the update ith iterationi=Ii-1+αiriSpecially：
The multi -components inverse migration of the update ith iteration records di=di-1+αiδdiSpecially：
Step S1200：After ith iteration is complete, the cost functional value f of ith iteration is calculatedi, whether judge current iteration
Meet convergence, it is final migrated section that newest migrated section is exported if meetingOtherwise it repeats
Step S700 to S1100, until obtaining final migrated section；
The cost functional value f for calculating ith iterationiSpecially：
In the equation (11), summation variable k includes tri- directions x, y and z of cartesian coordinate；
The convergence is specially：
In the equation (12), theshold indicates the threshold value standard of iteration stopping, usually chooses 0.001.
The present invention relates to a kind of high efficiency multi-component seismic data true amplitude migration imaging methods, and the present invention is in elastic wave
Improved on the basis of reverse-time migration method, be characterized in directly with multi-component seismic data be input, by encryption algorithm,
The measures such as segmentation imaging, then under the frame of inverting, realize the elastic wave reverse-time migration based on inverting.The present invention will encode, when
Between be segmented and the thought of inverting is introduced into elastic wave reverse-time migration, compared with common elastic wave reverse-time migration, the present invention can obtain
The pre-stack depth migration section of high-precision and hi-fi of amplitude is obtained, while the efficiency of algorithm is very high.The present invention is completely kept in length and breadth
Vectorial property, amplitude and the phase characteristic of wave, migrated section have fidelity；Crosstalk is made between effectively eliminating longitudinal and shear wave
At offset illusion, improve the precision of imaging；The amount of storage for effectively reducing elastic wave reverse-time migration improves calculating effect
Rate.
For the feasibility and validity further illustrated the present invention, three examples are named：
Example 1：
Fig. 2 is the hollow illustraton of model of two dimension；Wherein, Fig. 2 (a) is velocity of longitudinal wave model；Fig. 2 (b) is shear wave velocity model；Figure
2 (c) is background velocity of longitudinal wave model；Fig. 2 (d) is background shear wave velocity model.21 ray parameters are set on this model, are risen
Beginning ray parameter p=-117 μ s/m, p=15 μ s/m are divided between ray parameter, and focus used is hypocenter of the explosion, source wavelet setting
For Ricker wavelet, dominant frequency is 20 hertz.Fig. 3 is the coding big gun collection cut off after direct wave, corresponding ray parameter p=0 μ s/m：Its
In, Fig. 3 (a) is the observational record horizontal component of coding；Fig. 3 (b) is the observational record vertical component of coding；Fig. 3 (c) is to utilize
The initial inverse migration recording level component of coding obtained by the method for the present invention；Fig. 3 (d) is to utilize the coding obtained by the method for the present invention
Initial inverse migration record vertical component；Fig. 3 (e) is the final inverse migration recording level using the coding obtained by the method for the present invention
Component；Fig. 3 (f) is the final inverse migration record vertical component using the coding obtained by the method for the present invention.It can from Fig. 3
Go out, inverse migration record obtained by the method for the present invention matches very good with observational record, this demonstrates the validity of the method for the present invention.Figure
4 be more big guns superposition migrated section of the hollow model of two dimension shown in Fig. 2：Wherein, Fig. 4 (a) is to utilize the I obtained by conventional methodαIt cuts open
Face, Fig. 4 (b) are to utilize the I obtained by conventional methodβSection, Fig. 4 (c) are to utilize the I obtained by the method for the present inventionαSection, Fig. 4 (d)
It is to utilize the I obtained by the method for the present inventionβSection.As can be seen that two sections are in the presence of brighter from Fig. 4 (a) and Fig. 4 (b)
Aobvious crosstalk noise and illusion, the resolution ratio of section is relatively low, and amplitude is unbalanced, as can be seen that originally from Fig. 4 (c) and Fig. 4 (d)
Section obtained by inventive method has suppressed crosstalk noise and illusion well, migrated section have very high precision, resolution ratio and
Signal-to-noise ratio, and amplitude is fidelity.By Fig. 3 and Fig. 4 it can be proved that the offset that the method for the present invention can obtain high quality is cutd open
Face, while also demonstrating the feasibility and validity of the method for the present invention.
Example 2：
Fig. 5 is Marmousi-ii models：Wherein, Fig. 5 (a) is velocity of longitudinal wave model；Fig. 5 (b) is shear wave velocity model；
Fig. 5 (c) is background velocity of longitudinal wave model；Fig. 5 (d) is background shear wave velocity model.The model be the various offset methods of verification at
As one of the international standard model of effect.31 ray parameters of setting on this model, start ray parameter p=-333 μ s/m,
P=20 μ s/m are divided between ray parameter, focus used is hypocenter of the explosion, and source wavelet is set as Ricker wavelet, and dominant frequency is 15 hertz
Hereby.Fig. 6 is more big guns superposition migrated section of Marmousi-ii models shown in Fig. 5：Wherein, Fig. 6 (a) is to utilize conventional method institute
The I obtainedαSection；Fig. 6 (b) is to utilize the I obtained by conventional methodβSection；Fig. 6 (c) is to utilize the I obtained by the method for the present inventionαIt cuts open
Face；Fig. 6 (d) is to utilize the I obtained by the method for the present inventionβSection.From Fig. 6 (a) and Fig. 6 (b) as can be seen that obtained by conventional method
IαSection and IβSection scarce capacity in terms of suppressing longitudinal and shear wave crosstalk, and imaging section amplitude serious unbalance, deep section
Relative amplitude and infidelity.In contrast, imaging section longitudinal and shear wave crosstalk obtained by the method for the present invention is almost completely eliminated, and is deviated
The resolution ratio and precision higher of section, signal-to-noise ratio is also more preferable, amplitude equalization more preferably, therefore it may be concluded that side of the present invention
Method imaging effect is more preferable.
Example 3：
Fig. 7 is SEG/EAGE Salt models：Wherein, Fig. 7 (a) is velocity of longitudinal wave model；Fig. 7 (b) is shear wave velocity mould
Type；Fig. 7 (c) is background velocity of longitudinal wave model；Fig. 7 (d) is background shear wave velocity model.The model is the various offset methods of verification
One of international standard In A Salt-dome Model of imaging effect.21 ray parameters, start ray parameter p=-333 are set on this model
μ s/m are divided into p=30 μ s/m between ray parameter, and focus used is hypocenter of the explosion, and source wavelet is set as Ricker wavelet, and dominant frequency is
15 hertz.Fig. 8 is more big guns superposition migrated section of SEG/EAGESalt models shown in Fig. 7：Wherein, Fig. 8 (a) is to utilize tradition side
I obtained by methodαSection；Fig. 8 (b) is to utilize the I obtained by conventional methodβSection；Fig. 8 (c) is to utilize the I obtained by the method for the present inventionα
Section；Fig. 8 (d) is to utilize the I obtained by the method for the present inventionβSection.As can be seen from Figure 8, offset obtained by the method for the present invention is cutd open
Face has higher resolution ratio, better signal-to-noise ratio, precision also higher, we can also clearly find out the method for the present invention in addition
Gained section has very high hi-fi of amplitude.
Above description and examples are merely to illustrate the present invention, it is every carry out based on the technical solution of the present invention etc.
Same changes and improvements, should all be included in the protection scope of the present invention.
Claims (6)
1. a kind of high efficiency multi-component seismic data true amplitude migration imaging method, it is characterised in that include the following steps：
(a)：The more big gun multi component seismic records D (x obtained by field inspectionr,t；xs)=(Dx,Dy,Dz), give background mould
Type, parameter preset determine the multi component seismic records for carrying out migration imaging, observation system parameter and source function S (x, t；xs),
Wherein, x representation spaces position vector, xsIndicate hypocentral location vector, xrIndicate that geophone station position vector, t indicate time variable；
(b)：Observation system parameter based on reading, design plane wave coding function C (xs, p), specially：
C(xs, p) and=fδ(t-p·xs), (1)
Wherein, p is ray parameter variable, fδIndicate impulse function；Plane wave coding function C (x hereins, p) and it is equivalent to time shift letter
Number, time shift amount pxsWith focal point position xsDifferent and change, ray parameter p is specially：
Wherein, α is earth's surface velocity of longitudinal wave, and θ is incidence angle；By different incidence angle θs, different ray parameter p can be obtained；
By the coding function, and using each ray parameter p to more big gun multi component seismic records D (x of observationr,t；xs)
It is encoded, obtains the coding multi -components big gun collection D about each ray parameter pC(xr,t；P), specially：
DC(xr,t；P)=∫ D (xr,t；xs)·C(xs,p)dxs； (3)
Equally, using identical coding function, and using each ray parameter p in observation system described in step (a)
Source function S (x, t；xs) encoded, obtain the coded source function S of each ray parameter pC(x,t；P), specially：
SC(x,t；P)=∫ S (x, t；xs)·C(xs,p)dxs； (4)
(c)：Utilize background model and the dominant record time t of more big gun multi component seismic records of observationmax, design record time t
Piecewise interval, specially：It will record time t point be three sections, be indicated with following formula：
Wherein, tR1, tR2And tR3Three interval ranges of record time, t are indicated respectively1(x) and t2(x) specific calculating formula is：
Wherein, αmaxFor the maximum speed of background model, hzRepresent depth, hxRepresent horizontal distance；
(d)：For each ray parameter p, it is performed both by following operation：Using each ray parameter p, by step (b)
Coding mode obtain corresponding with ray parameter coded source function, the elastic wave side of the isotropic medium based on routine
It is corresponding every to obtain the ray parameter into traveling-wave field forward direction continuation to the coded source using background model given herein above for journey
The source wavefield at a moment utilizes the record time slice in step (c) in each time step of the source wavefield continuation
Section is as time constraint condition, using conventional elastic wave field energy balane equation calculation source wavefield each moment, each
The wavefield energy value of position, is judged by energy values size, is obtained in each time interval, each position wavefield energy maximum
At the time of value is corresponding, using this moment as the imaging time T (x) in the time interval, corresponding to the position, the moment,
At the position, using image-forming condition formula, preserve in image-forming condition with the relevant item of source wavefield, while preserve the moment, should
The wave field maximum energy value of position；
And then it obtains in the corresponding imaging time of each ray parameter, image-forming condition with the relevant item of source wavefield and wave field most
Big energy value；Finally the wave field maximum energy value in each corresponding period of each ray parameter is overlapped, is obtained
Obtain illuminance compensation factor；
(e)：For each ray parameter p, it is performed both by following operation：To each ray parameter p, equally based on described each
To the equations for elastic waves of same sex medium, to step (b) the coding multi -components big gun collection D that ray parameter p is correspondingC(xr,t；p)
Into the continuation of traveling-wave field inverse time, the geophone station wave field at each corresponding moment of the ray parameter is obtained, is prolonged in the geophone station wave field
In each time step opened up, using the imaging time T (x) of each period, each position that are preserved in step (d), judge
Whether current time meets imaging time, if conditions are not met, then without imaging, if it is satisfied, then according to image-forming condition, utilizes
It is imaged with the relevant item of source wavefield and geophone station wave field in the image-forming condition preserved in step (d), obtains the ray
The corresponding migrated section of parameter；
And then obtain the corresponding migrated section of each ray parameter；All corresponding migrated sections of ray parameter are folded
Add, while energy compensating is carried out to the migrated section of superposition using illuminance compensation factor in step (d), obtains initial offset and cut open
The migrated section of face namely the 0th iteration
(f)：For each ray parameter p, following operation is executed：To each ray parameter p, pass through the volume in step (b)
Code mode obtains coded source function corresponding with the ray parameter, equally the elastic wave side based on the isotropic medium
Journey obtains the corresponding back of the body of the ray parameter to the coded source using background model given herein above into traveling-wave field forward direction continuation
Scape wave field；Utilize the initial offset section I described in step (e)0, the ray parameter background wave field and step
(a) background model given in, builds the corresponding virtual focus of the ray parameter；Based on each to same of the virtual focus
Property dielectric resilient wave equation, using given background model to the corresponding coded source function of the ray parameter into traveling-wave field forward direction
Continuation obtains the corresponding inverse migration wave field of the ray parameter, and record sampling is carried out to the corresponding inverse migration wave field of the ray parameter,
Obtain the corresponding coding multi -components inverse migration big gun collection of the ray parameter；
And then obtain the corresponding coding multi -components inverse migration big gun collection of each ray parameter, the coding multi -components big gun collection structure of all big guns
At initial code multi -components inverse migration big gun collection namely the coding multi -components inverse migration big gun collection of the 0th iterationAt this point, it is i=0 that current iterations, which are arranged,；
(g)：Current iterations i=i+1 is set, based on the migrated section obtained by (i-1)-th iterationWith the coding multi -components inverse migration big gun collection obtained by (i-1)-th iterationCarry out inverting iteration；
In ith iteration, for each ray parameter p, by the corresponding coding of the ray parameter obtained by (i-1)-th iteration
Multi -components inverse migration big gun collection di-1The coding multi -components big gun collection D of observation corresponding with the ray parameterCIt is poor to make, and calculates ray ginseng
The corresponding coding multi -components big gun collection residual error of number；And then obtain the corresponding coding multi -components big gun collection residual error of all ray parameters
It is described by the corresponding coding multi -components inverse migration big gun collection d of the ray parameter obtained by (i-1)-th iterationi-1With the ray
The coding multi -components big gun collection D of the corresponding observation of parameterCIt is poor to make, and calculates the corresponding coding multi -components big gun collection residual error of the ray parameter,
Specially：
Δdi-1=di-1-DC； (7)
(h)：For each ray parameter p, following operation is executed：To each ray parameter p, equally based on it is described it is each to
The equations for elastic waves of same sex medium, coding multi -components big gun collection residual error corresponding to the ray parameterInto the continuation of traveling-wave field inverse time, each corresponding moment of the ray parameter is obtained
Geophone station wave field, in each time step of the geophone station wave field extrapolation, using each period preserved in step (d),
The imaging time T (x) of each position, judges whether current time meets imaging time, if conditions are not met, then without being imaged,
If it is satisfied, then according to image-forming condition, using in the image-forming condition preserved in step (d) with the relevant item of source wavefield and inspection
Wave point wave field is imaged, and the corresponding gradient profile of the ray parameter is obtained；
And then obtain the corresponding gradient profile of each ray parameter；All corresponding gradient profiles of ray parameter are folded
Add, while energy compensating is carried out to the gradient profile of superposition using illuminance compensation factor in step (d), obtains current iteration
The gradient profile of gradient profile namely ith iteration
(i)：Using optimal inversion algorithm, the descent direction section of ith iteration is built
(j)：For each ray parameter p, coding corresponding with the ray parameter is obtained by the coding mode in step (b)
Source function, the equally equations for elastic waves based on the isotropic medium, using background model given herein above to the volume
Code focus obtains the corresponding background wave field of the ray parameter into traveling-wave field forward direction continuation；Utilize i-th obtained in step (i)
The descent direction section r of secondary iterationi, the ray parameter background wave field and step (a) in give background model,
Build the corresponding virtual focus of the ray parameter；Equally based on the isotropic medium elastic wave side with the virtual focus
Journey obtains this and penetrates to the corresponding coded source function of the ray parameter using given background model into traveling-wave field forward direction continuation
The corresponding disturbance inverse migration wave field of line parameter, disturbance inverse migration wave field corresponding to the ray parameter carry out record sampling, obtain
The corresponding coding multi -components of the ray parameter disturb inverse migration big gun collection；
And then the corresponding coding multi -components disturbance inverse migration big gun collection of each ray parameter is obtained, the coding multi -components of all big guns are disturbed
Dynamic big gun collection constitutes the coding multi -components disturbance inverse migration big gun collection of ith iterationIt is sharp again
The optimization step-length α of current iteration is calculated with stepsize formulai；
The stepsize formula is
In the equation (8), summation variable k includes tri- directions x, y and z of cartesian coordinate,Indicate ith iteration
Encode multi -components disturbance inverse migration big gun collection δ diK durection components,Indicate the coding multi -components big gun collection residual delta of ith iteration
diK durection components, k=x, y, z；
(k)：The optimization step-length α obtained using step (j)iAnd the descent direction section r that step (i) obtainsi, update ith iteration
Migrated section Ii=Ii-1+αiri；Optimize step-length α using statingiAnd the coding multi -components disturbance inverse migration big gun collection that step (j) obtains
δdi, update the coding multi -components inverse migration big gun collection d of ith iterationi=di-1+αiδdi；
The migrated section I of the update ith iterationi=Ii-1+αiriSpecially：
The multi -components inverse migration of the update ith iteration records di=di-1+αiδdiSpecially：
(l)：After ith iteration is complete, the cost functional value f of ith iteration is calculatedi, judge whether current iteration meets convergence mark
Standard, it is final migrated section that newest migrated section is exported if meetingOtherwise step (g) is repeated extremely
(k), until obtaining final migrated section；
The cost functional value f for calculating ith iterationiSpecially：
In the equation (11), summation variable k includes tri- directions x, y and z of cartesian coordinate；
The convergence is specially：
In the equation (12), theshold indicates the threshold value standard of iteration stopping, usually chooses 0.001.
2. a kind of high efficiency multi-component seismic data true amplitude migration imaging method according to claim 1, feature exist
It is in the equations for elastic waves of step (d), step (e), step (f), the isotropic medium described in step (h) and step (j)：
In the equation (13), V=(Vx,Vy,Vz)TIndicate vector velocity wave field；σ=(σxx,σyy,σzz,σxy,σyz,σxz)TIt indicates
Vector stress wave field；λ and μ indicates medium Lame constants；ρ indicates Media density；X, y and z indicate the x of cartesian coordinate, y respectively
With tri- directions z.
3. a kind of high efficiency multi-component seismic data true amplitude migration imaging method according to claim 2, feature exist
In step (d), the image-forming condition described in step (e) and step (h) is：
Wherein, α is medium velocity of longitudinal wave, and β is medium shear wave velocity, the p-and s-wave velocity of medium and the medium Lame constants,
Meet relationship between Media density：
With
In the equation (14), in step (e),In step (h),
In the equation (14),Indicate the vector velocity wave field in source wavefield,Indicate the vector stress wave field in geophone station wave field.
4. a kind of high efficiency multi-component seismic data true amplitude migration imaging method according to claim 3, feature exist
Described in step (d) in each time step of the source wavefield continuation, utilize the record time slice area in step (c)
Between be used as time constraint condition, using elastic wave field energy balane equation calculation each moment, the wavefield energy value of each position,
Specially：Elastic wave field energy balane equation is：
E=K+W, (15)
Wherein, E is the energy of elastic wave field, and K is the kinetic energy of elastic wave field, and W is the potential energy of elastic wave field；The kinetic energy of elastic wave field
It can be obtained by its velocity field and stress field calculation with potential energy, specially：
Wherein, U is displacement field, which can be obtained by doing time integral to the vector velocity wave field V；
Judged by energy values size described in step (d), obtain in each time interval, each position wavefield energy most
At the time of corresponding to big value, using this moment as the imaging time T (x) in the time interval, corresponding to the position, specially：
With first time interval t in step (c)R1For, in each position, if current time, t was in tR1Within, then when
The wavefield energy E at preceding momentt(x) if more than in the time interval, the maximum value E of the wavefield energy at all moment beforemax(x),
So, Emax(x)=Et(x), the imaging time T (x) of the time interval=t；Otherwise, in the time interval, the wave field of the position
The maximum value and imaging time of energy remain unchanged；For other two time interval, then similarly；
Described in step (d) at the moment, the position, using image-forming condition formula, preserve in image-forming condition with focus wave
Relevant item, at the same preserve the moment, the position wave field maximum energy value, specially：In image-forming condition with source wavefield
Relevant item includes：WithIt totally nine, calculates
This nine value and the value of calculating is preserved at the moment, the position；Meanwhile preserve the moment, the position wave field most
Big energy value Emax(x)。
5. a kind of high efficiency multi-component seismic data true amplitude migration imaging method according to claim 3, feature exist
The imaging time T of each period, each position that are preserved in utilization step (d) described in step (e) and step (h)
(x), judge whether current time meets imaging time, if conditions are not met, then without imaging, if it is satisfied, then according to imaging
Condition is imaged with the relevant item of source wavefield and geophone station wave field using in the image-forming condition preserved in step (d), is obtained
The corresponding migrated section of the ray parameter is obtained, specially：Meet first time interval t in step (c) with current time tR1For
Example, in each position, if current time t satisfaction：| T (x)-t | < 5, then according to the image-forming condition, to the focus of preservation
The continuous item of geophone station wave field obtained by wave field continuous item and current time applies the image-forming condition, obtains m1And m2, and should
The imaging results at moment are added in the imaging section of the position；Otherwise, the moment, the position are without imaging；For in addition
Two time intervals, then similarly.
6. a kind of high efficiency multi-component seismic data true amplitude migration imaging method according to claim 5, feature exist
Based on the isotropic medium equations for elastic waves with the virtual focus described in step (f) and step (j), using given
Background model to the corresponding coded source function of the ray parameter into traveling-wave field forward direction continuation, herein carry the virtual shake
The isotropic medium equations for elastic waves in source is expressed as：
In the equation (17), δ V=(δ Vx,δVy,δVz)TIndicate inverse migration vector velocity wave field；δ σ=(δ σxx,δσyy,δσzz,δ
σxy,δσyz,δσzx)TIndicate inverse migration vector stress wave field；
In the equation (17), virtual source vector F=(Fxx,Fyy,Fzz,Fxy,Fyz,Fxz)TIt is expressed as：
In the equation (18), intermediate variable m '1With m '2With the m1And m2Meet following relationship：
With
In step (f),In step (j),
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201611197758.4A CN106842300B (en) | 2016-12-21 | 2016-12-21 | A kind of high efficiency multi-component seismic data true amplitude migration imaging method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201611197758.4A CN106842300B (en) | 2016-12-21 | 2016-12-21 | A kind of high efficiency multi-component seismic data true amplitude migration imaging method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106842300A CN106842300A (en) | 2017-06-13 |
| CN106842300B true CN106842300B (en) | 2018-10-30 |
Family
ID=59136914
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201611197758.4A Expired - Fee Related CN106842300B (en) | 2016-12-21 | 2016-12-21 | A kind of high efficiency multi-component seismic data true amplitude migration imaging method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106842300B (en) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108241173B (en) * | 2017-12-28 | 2019-06-21 | 中国石油大学(华东) | A kind of seismic data offset imaging method and system |
| CN111999770B (en) * | 2020-09-03 | 2024-01-16 | 中国地质科学院 | TTI medium conversion PS wave precise beam offset imaging method and system |
| CN112462428B (en) * | 2020-11-13 | 2024-02-20 | 中国地质科学院 | Multi-component seismic data migration imaging method and system |
| CN112904426B (en) * | 2021-03-27 | 2022-09-30 | 中国石油大学(华东) | Decoupling elastic wave reverse time migration method, system and application |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102116870A (en) * | 2011-02-12 | 2011-07-06 | 中国石油大学(华东) | Elastic wave gaussian beam pre-stack depth migration technology |
| CN102156296A (en) * | 2011-04-19 | 2011-08-17 | 中国石油大学(华东) | Elastic reverse time migration imaging method by combining seismic multi-component |
| CN103760603A (en) * | 2014-01-28 | 2014-04-30 | 中国石油大学(北京) | Pre-stack time migration method and device for converted wave seismic data |
| CN104991268A (en) * | 2015-07-03 | 2015-10-21 | 中国地质大学（北京） | True amplitude migration imaging method |
| CN105388520A (en) * | 2015-10-22 | 2016-03-09 | 中国石油化工股份有限公司 | Seismic data pre-stack reverse time migration imaging method |
| CN105807315A (en) * | 2016-03-14 | 2016-07-27 | 中国石油大学(华东) | Elastic vector reverse time migration imaging method |
| CN105974470A (en) * | 2016-07-04 | 2016-09-28 | 中国石油大学(华东) | Multi-component seismic data least squares reverse time migration imaging method and system |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6819628B2 (en) * | 2003-04-07 | 2004-11-16 | Paradigm Geophysical (Luxembourg) S.A.R.L. | Wave migration by a krylov space expansion of the square root exponent operator, for use in seismic imaging |
| CN101285894B (en) * | 2008-05-30 | 2011-02-09 | 中国科学院地质与地球物理研究所 | Heaved earth surface collected seismic data direct prestack time migration method |
| CN101957455B (en) * | 2010-09-20 | 2012-07-18 | 中国海洋石油总公司 | Method of three-dimensional preserved-amplitude pre-stack time migration |
| CN103207409B (en) * | 2013-04-17 | 2016-01-06 | 中国海洋石油总公司 | A kind of frequency domain full-waveform inversion seismic velocity modeling method |
-
2016
- 2016-12-21 CN CN201611197758.4A patent/CN106842300B/en not_active Expired - Fee Related
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102116870A (en) * | 2011-02-12 | 2011-07-06 | 中国石油大学(华东) | Elastic wave gaussian beam pre-stack depth migration technology |
| CN102156296A (en) * | 2011-04-19 | 2011-08-17 | 中国石油大学(华东) | Elastic reverse time migration imaging method by combining seismic multi-component |
| CN103760603A (en) * | 2014-01-28 | 2014-04-30 | 中国石油大学(北京) | Pre-stack time migration method and device for converted wave seismic data |
| CN104991268A (en) * | 2015-07-03 | 2015-10-21 | 中国地质大学（北京） | True amplitude migration imaging method |
| CN105388520A (en) * | 2015-10-22 | 2016-03-09 | 中国石油化工股份有限公司 | Seismic data pre-stack reverse time migration imaging method |
| CN105807315A (en) * | 2016-03-14 | 2016-07-27 | 中国石油大学(华东) | Elastic vector reverse time migration imaging method |
| CN105974470A (en) * | 2016-07-04 | 2016-09-28 | 中国石油大学(华东) | Multi-component seismic data least squares reverse time migration imaging method and system |
Non-Patent Citations (3)
| Title |
|---|
| "基于时移成像条件的波动方程保幅成像";李振春 等;《勘探地球物理进展》;20080430;第31卷(第4期);第270-273页 * |
| "平面波最小二乘逆时偏移编码策略分析";李闯 等;《石油物探》;20150930;第54卷(第5期);第592-601页 * |
| "真振幅偏移方法综述";李振春 等;《勘探地球物理进展》;20080229;第31卷(第1期);第10-15页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN106842300A (en) | 2017-06-13 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105974470B (en) | A kind of multi-component seismic data least square reverse-time migration imaging method and system | |
| CN108873066B (en) | Elastic medium wave equation reflected wave travel time inversion method | |
| CN106842300B (en) | A kind of high efficiency multi-component seismic data true amplitude migration imaging method | |
| US10234582B2 (en) | Joint inversion of seismic data | |
| CN108241173B (en) | A kind of seismic data offset imaging method and system | |
| EP2946233B1 (en) | Method for ray based tomography guided by waveform inversion | |
| CN109557582B (en) | A kind of two dimension multi-component seismic data offset imaging method and system | |
| US11092708B2 (en) | Processes and systems to enhance illumination and resolution of seismic images using multiple reflected wavefields | |
| CN105652321A (en) | Visco-acoustic anisotropic least square inverse time migration imaging method | |
| CN109856679B (en) | Method and system for imaging elastic wave Gaussian beam offset of anisotropic medium | |
| CN110780351B (en) | Longitudinal wave and converted wave prestack joint inversion method and system | |
| CN111045077B (en) | Full waveform inversion method of land seismic data | |
| RU2570827C2 (en) | Hybrid method for full-waveform inversion using simultaneous and sequential source method | |
| US11635540B2 (en) | Methods and devices performing adaptive quadratic Wasserstein full-waveform inversion | |
| Ibrahim | Separating simultaneous seismic sources using robust inversion of Radon and migration operators | |
| CN109143352B (en) | A kind of anisotropic medium Seismic reflection character establishing equation method | |
| Bai et al. | Turning high resolution FWI model into pseudo-reflectivity: A case study using sparse OBN | |
| US11604299B2 (en) | Mixed-phase source wavelet estimation from recorded seismic data | |
| Guo et al. | Becoming effective velocity-model builders and depth imagers, Part 2—The basics of velocity-model building, examples and discussions | |
| Jin et al. | 2D multiscale non‐linear velocity inversion | |
| Vinje et al. | Preserved‐traveltime smoothing | |
| CN108375794A (en) | Based on the VSP fracture hole Diffraction Imaging technical methods symmetrically observed | |
| CN109901221B (en) | Seismic data anisotropy modeling method based on dynamic correction velocity parameter | |
| Wu et al. | Attenuation compensation in multicomponent Gaussian beam prestack depth migration | |
| Kazei et al. | Acquisition and near-surface impacts on VSP mini-batch FWI and RTM imaging in desert environment |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20181030 |