CN108363097A - A kind of seismic data offset imaging method - Google Patents
A kind of seismic data offset imaging method Download PDFInfo
- Publication number
- CN108363097A CN108363097A CN201810104681.4A CN201810104681A CN108363097A CN 108363097 A CN108363097 A CN 108363097A CN 201810104681 A CN201810104681 A CN 201810104681A CN 108363097 A CN108363097 A CN 108363097A
- Authority
- CN
- China
- Prior art keywords
- domain
- vertical time
- depth
- time
- coordinate
- 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.)
- Granted
Links
- 238000003384 imaging method Methods 0.000 title claims abstract description 38
- 238000013508 migration Methods 0.000 claims abstract description 69
- 230000005012 migration Effects 0.000 claims abstract description 69
- 238000000034 method Methods 0.000 claims abstract description 39
- 238000005457 optimization Methods 0.000 claims abstract description 17
- 238000006243 chemical reaction Methods 0.000 claims description 17
- 230000009466 transformation Effects 0.000 claims description 11
- 239000002245 particle Substances 0.000 claims description 6
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 3
- 241001269238 Data Species 0.000 claims description 3
- 238000012952 Resampling Methods 0.000 claims description 3
- 238000013213 extrapolation Methods 0.000 claims description 3
- 238000005070 sampling Methods 0.000 claims description 3
- 238000006073 displacement reaction Methods 0.000 claims 1
- 238000007689 inspection Methods 0.000 claims 1
- 238000004364 calculation method Methods 0.000 abstract description 4
- 238000005516 engineering process Methods 0.000 description 8
- 238000007796 conventional method Methods 0.000 description 5
- 238000010586 diagram Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 230000008901 benefit Effects 0.000 description 2
- 241000208340 Araliaceae Species 0.000 description 1
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 1
- 235000003140 Panax quinquefolius Nutrition 0.000 description 1
- 238000004880 explosion Methods 0.000 description 1
- 235000008434 ginseng Nutrition 0.000 description 1
- 230000000750 progressive effect Effects 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
- G01V1/282—Application of seismic models, synthetic seismograms
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
The present invention discloses a kind of seismic data offset imaging method.This method includes：Obtain seismic data, observation system parameter, Migration velocity model and the offset parameter of pending migration imaging；Migration velocity model is transformed into vertical time-domain by Depth Domain；Isotropic medium one-order velocity stress ACOUSTIC WAVE EQUATION is transformed into vertical time-domain by Depth Domain；Obtain the speed stress optimization operator per each corresponding moment of big gun；Source wavefield reconstruction is carried out using the speed stress optimization operator at each corresponding moment of the big gun of gained, obtains the source wavefield at each moment of the big gun；Obtain the geophone station wave field at each moment per big gun；It obtains per the corresponding single-shot migrated section of big gun；The single-shot migrated section of all big guns is overlapped, final migrated section is obtained.The method of the present invention can effectively reduce the calculation amount and amount of storage of conventional three-dimensional reverse-time migration method, improve and calculate efficiency.
Description
Technical field
The invention belongs to seismic exploration technique fields, more particularly to a kind of seismic data offset imaging method.
Background technology
Seismic data reverse-time migration imaging technique is based on round trip wave theory, is not limited by inclination angle and lateral velocity variation, gram
The limitation for having taken traditional radiographic class offset and one-way wave offset method has good in complex oil and gas reservoir imaging side face
Application prospect and economic benefit.
For reverse-time migration method, the storage of source wavefield is the bottleneck for restricting this method and applying in practice.In order to gram
The storage problem is taken, common several method technology is in the industry at present：Checkpoint technologies, RANDOM BOUNDARY wave-field reconstruction skill
Art, efficiency frontier memory technology.These types of technology has respective limitation：Checkpoint technology calculation amounts are larger；At random
Boundary wave-field reconstruction technique, which introduces additional random noise, influences image quality；Efficiency frontier memory technology is for three dimensions
According to amount of storage is still difficult to receive.
For the amount of storage and computational efficiency problem present in current reverse-time migration method, it is necessary to establish it is a set of it is new, with
Seismic data offset imaging method with the characteristics of " high efficiency, high-precision, low storage ".
Invention content
The purpose of the present invention is to provide a kind of seismic data offset imaging methods, to improve seismic data migration imaging
Efficiency.
To achieve the above object, the present invention provides following schemes：
A kind of seismic data offset imaging method, the method includes：
Earthquake work area observation system parameter, rate pattern and more big gun seismic datas of observation are read, determines and carries out partially
Move into more big guns observation earthquake record, observation system parameter, Depth Domain Migration velocity model and the offset parameter of picture；
Using the vertical time-domain coordinate conversion relation formula of Depth Domain-, Depth Domain coordinate is transformed to vertical time-domain, is based on
The Depth Domain Migration velocity model obtains the Migration velocity model of vertical time-domain, in turn using cubic spline interpolation algorithm
Obtain the smoothed offset rate pattern of vertical time-domain；
Using the vertical time-domain coordinate conversion relation formula of Depth Domain-, the isotropic medium one-order velocity-of Depth Domain is answered
Power ACOUSTIC WAVE EQUATION is transformed to isotropic medium one-order velocity-stress ACOUSTIC WAVE EQUATION of vertical time-domain；
Shot position coordinate according to the observation system parameter acquiring per big gun；It is set respectively at each shot position coordinate
Set source wavelet；According to the vertical time-domain Migration velocity model, vertical time-domain smoothed offset rate pattern and offset ginseng
Number solves the vertical time-domain single order stress-speed ACOUSTIC WAVE EQUATION using the preferential difference numerical algorithm of staggered-mesh, realizes
Per the positive continuation of big gun source wavefield, the source wavefield at each moment per big gun is obtained；Preserve the speed at each moment per big gun
Degree-stress optimization operator；
Geophone station position coordinates according to the observation system parameter acquiring per big gun；According to the vertical time-domain offset speed
Model, vertical time-domain smoothed offset rate pattern and offset parameter are spent, the vertical time-domain is solved using numerical algorithm
Single order stress-speed ACOUSTIC WAVE EQUATION carries out inverse time continuation to the observation earthquake record of every big gun；When obtaining each per big gun
The geophone station wave field at quarter；
Read storage every big gun each moment speed-stress optimization operator, using reading Optimizing operator and
Specific difference scheme carries out source wavefield reconstruction, obtains the source wavefield of the reconstruction at each moment per big gun；
Mutually in the same time, it is being based on migration imaging condition, to the source wavefield of the reconstruction at each moment of every big gun
It is imaged with geophone station wave field, obtains the single-shot migrated section per big gun；And then the single-shot migrated section of all big guns is obtained, by institute
There is the single-shot migrated section of big gun to be overlapped, obtains the final migrated section of vertical time-domain；
Using the vertical time-domain coordinate inverse transformation relationship of Depth Domain-, by the final migration imaging of the vertical time-domain
Section is converted into Depth Domain by vertical time-domain, obtains the final migrated section of Depth Domain.
Optionally, described to utilize the vertical time-domain coordinate conversion relation formula of Depth Domain-, Depth Domain coordinate is transformed to vertical
Time-domain is based on the Depth Domain Migration velocity model, using cubic spline interpolation algorithm, obtains the offset speed of vertical time-domain
Model is spent, and then obtains the smoothed offset rate pattern of vertical time-domain, specially：
The vertical time-domain coordinate conversion relation formula of Depth Domain-is：
In formula (1), x, y and z indicate the x, y and z directionss of cartesian coordinate system respectively；τ (x, y, z) is by Depth Domain coordinate
The vertical time-domain coordinate variable of variable z conversion gained, it is one-to-one with z；It is by the Depth Domain
Migration velocity model vdThe Depth Domain smoothed offset rate pattern of (x, y, z) smooth gained；Depth Domain can be sat using formula (1)
Mark is transformed to vertical time-domain；
It is described that vertical time-domain is obtained using cubic spline interpolation algorithm based on the Depth Domain Migration velocity model
Migration velocity model, and then the smoothed offset rate pattern of vertical time-domain is obtained, it specifically includes：
(a) according to the one-to-one relationship between Depth Domain coordinate z and vertical time-domain coordinate τ, the Depth Domain is smooth
Migration velocity modelThe Migration velocity model v of as original vertical time-domainτ1(x, y, τ), at this time vertical time
The coordinate interval of domain coordinate τ is non-uniform；
(b) the maximum coordinates interval for the vertical time-domain being subject to described in (a), using cubic spline interpolation algorithm to institute
The vertical time-domain coordinate stated is into row interpolation resampling, when obtaining uniform vertical time-domain coordinate, while can be obtained vertical
Between domain Migration velocity model vτ(x,y,τ)；
Similarly, by Depth Domain smoothed offset rate patternThe smoothed offset speed of vertical time-domain can be obtained
Spend model
Optionally, described to utilize the vertical time-domain coordinate conversion relation formula of Depth Domain-, by the isotropic medium of Depth Domain
One-order velocity-stress ACOUSTIC WAVE EQUATION is transformed to isotropic medium one-order velocity-stress ACOUSTIC WAVE EQUATION of vertical time-domain, specifically
For：
Depth Domain isotropic medium one-order velocity-stress ACOUSTIC WAVE EQUATION is specially：
In formula (2), ux, uyAnd uzRespectively Depth Domain isotropic medium sound wave particle velocity wave field is in x, y and z tri-
The component in direction, p are Depth Domain isotropic medium sound wave direct stress wave field；
It is described by isotropic medium one-order velocity-stress ACOUSTIC WAVE EQUATION of Depth Domain be transformed to vertical time-domain it is each to
Same sex medium one-order velocity-stress ACOUSTIC WAVE EQUATION, specially：
(a) Depth Domain coordinate system (x, y, z) is transformed to vertical time-domain coordinate systemThe seat of two coordinate systems
Mark variable meets following relationship：
(b) created symbol fxIt indicates：
Define vertical time-domain operator α=τxWith β=τy, according to chain rule, the partial derivative of cartesian coordinate system with it is vertical
The partial derivative of time-domain coordinate system meets following relationship：
So unit base vector of cartesian coordinate systemWith the unit base vector of bent coordinate systemBetween
Relationship is：
(c) the arbitrary vector under cartesian coordinate systemIt is in vertical time-domain corresponding form
For：
Wherein, Vx, VyAnd VzThe component in tri- directions x, y and z is tied up in cartesian coordinate for vector V；Vξ,And VτFor arrow
V is measured in vertical time-domain coordinate system ξ,With the component in tri- directions τ；
With divergence theorem, by the component of its vertical time-domain of the gradient of the vector V under cartesian coordinate system and divergence
And unit base vector is indicated, specially：
(d) by formula (8), (9) substitute into formula (2), you can obtain the vertical time domain one-order velocity-stress of isotropic medium
ACOUSTIC WAVE EQUATION, specially：
In formula (10), Uξ,UτRespectively vertical time-domain isotropic medium sound wave particle velocity wave field is when vertical
Between domain coordinate system ξ,The component in tri- directions τ, P are vertical time-domain isotropic medium sound wave direct stress wave field.
Optionally, the speed-stress optimization operator for preserving each moment per big gun, specially：
Arbitrary t moment, speed-stress optimization operator AtIt is defined as follows：
In formula (11),It is just answered for the vertical time-domain sound wave on i-th of mesh point outside t moment model area
Reeb；Wherein i=1,2 ..., N indicate that half exponent number of space difference used in wave field extrapolation, N indicate maximum half exponent number；H is indicated
Space lattice spacing；aiFor first group of Optimizing operator coefficient.
Optionally, described to carry out source wavefield reconstruction using the Optimizing operator read and specific difference scheme, it obtains
The source wavefield of the reconstruction at each moment per big gun, specially；
By taking direct stress wave field P as an example, the space in vertical time-domain isotropic medium single order stress-speed ACOUSTIC WAVE EQUATION
The staggering mesh finite-difference form of partial derivative is：By taking P is in x=0 (i.e. left margin) as an example, i-th point on the right side of x=0 of P
First-order partial derivative It can be indicated with finite difference, specially：
Wherein, bi,jFor second group of Optimizing operator coefficient, the as described specific difference scheme of formula (12), by formula (12)
With formula (11) comparison it is found that formula (12) contains the Optimizing operator At, wave field can be found out on a left side by formula (12) formula
The partial derivative of N-1 points in boundary；
Similarly, the partial derivative of N-1 points in model its coboundary can be obtained；
For the point inside model, its partial derivative is calculated using conventional staggering mesh finite-difference, to realize focus
Wave-field reconstruction obtains the source wavefield of the reconstruction at each moment per big gun.
Optionally, described mutually in the same time, to be based on migration imaging condition, the reconstruction to each moment of every big gun
Source wavefield and geophone station wave field be imaged, obtain per big gun single-shot migrated section；And then the single-shot for obtaining all big guns is inclined
Section is moved, the single-shot migrated section of all big guns is overlapped, obtains the final migrated section of vertical time-domain, specially：
The image-forming condition is：
In formula in (13), PS(x, y, τ, t) and PR(x, y, τ, t) indicates the arbitrary spatial position at each moment respectively
Vertical time-domain source wavefield and geophone station wave field, Is(x, y, τ) indicates that the single-shot migrated section of s big guns, T indicate dominant record
Time；
It is described to be overlapped the single-shot migrated section of all big guns, the final migrated section of vertical time-domain is obtained, is had
Body is：
In formula (14), I (x, y, τ) indicates the final migrated section of vertical time-domain.
Optionally, described to utilize the vertical time-domain coordinate inverse transformation relationship of Depth Domain-, by the final of the vertical time-domain
Migration imaging section Depth Domain is converted by vertical time-domain, obtain the final migrated section of Depth Domain, specially：
The vertical time-domain coordinate inverse transformation relationship of Depth Domain-is：
In formula (15), Z (x, y, τ) is by the Depth Domain coordinate variable of vertical time-domain coordinate τ inverse transformation Depth Domains；At this time
Z (x, y, τ) samples for unequal interval, and equal interval sampling is obtained into row interpolation to Z (x, y, τ) using cubic spline interpolation
Depth Domain coordinate variable
Using formula (15), by the final migration imaging section I (x, y, τ) of the vertical time-domain by vertical time-domain
It is converted into Depth Domain, obtains the final migrated section of Depth Domain
The invention discloses technique effects once：
1) present invention is by introducing speed-stress optimization operator, and compared with conventional reverse-time migration method, the present invention effectively drops
The low amount of storage of reverse-time migration method, has pushed level of application of the reverse-time migration method in real data；2) present invention profit
With coordinate conversion relation, Depth Domain coordinate is converted into vertical time-domain coordinate, further reduced the meter of reverse-time migration method
Calculation amount improves the computational efficiency of reverse-time migration method；3) the method for the present invention is improving computational efficiency, is reducing the same of amount of storage
When, the precision of method has been effectively kept, the migration imaging section of high quality is can get；4) it present invention can be widely used to oil gas
It is more obvious especially for the imaging effect in complicated structure deep in Exploration Domain.
Description of the drawings
It in order to more clearly explain the embodiment of the invention or the technical proposal in the existing technology, below will be to institute in embodiment
Attached drawing to be used is needed to be briefly described, it should be apparent that, the accompanying drawings in the following description is only some implementations of the present invention
Example, for those of ordinary skill in the art, without having to pay creative labor, can also be according to these attached drawings
Obtain other attached drawings.
Fig. 1 is the flow diagram of seismic data offset imaging method of the present invention；
Fig. 2 is the hollow Migration velocity model figure of two dimension provided by the invention；Wherein, Fig. 2 (a) is depth domain model, Fig. 2
(b) it is vertical time domain model；
Fig. 3 is the single-shot earthquake record of the hollow model of two dimension shown in Fig. 2；Wherein, Fig. 3 (a) is the list obtained by conventional method
Big gun earthquake record, Fig. 3 (b) are the single-shot seismograms 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 obtained by conventional method
Depth Domain migrated section, Fig. 4 (b) is the vertical time-domain migrated section obtained by the method for the present invention, and Fig. 4 (c) is present invention side
Depth Domain migrated section obtained by method.
Specific implementation mode
With reference to the attached drawing in the embodiment of the present invention, technical solution in the embodiment of the present invention carries out clear, complete
Ground describes, it is clear that described embodiment is only a part of the embodiment of the present invention, instead of all the embodiments.It is based on
The embodiment of the present invention, those of ordinary skill in the art are obtained every other without making creative work
Embodiment shall fall within the protection scope of the present invention.
It is below in conjunction with the accompanying drawings and specific real to keep above-mentioned purpose, the features and advantages of the present invention more obvious and easy to understand
Applying mode, the present invention is described in further detail.
Fig. 1 is the flow diagram of seismic data offset imaging method of the present invention.As shown in Figure 1, the method includes：
Step 101：Earthquake work area observation system parameter, rate pattern and more big gun seismic datas of observation are read, really
Surely more big guns observation earthquake record, observation system parameter, Depth Domain Migration velocity model and the offset parameter of migration imaging are carried out.This
Step is initialization procedure, carries out subsequent processing on this basis.
Step 102：Migration velocity model is transformed into vertical time-domain by Depth Domain.Utilize the vertical time-domain of Depth Domain-
Depth Domain coordinate is transformed to vertical time-domain by coordinate conversion relation formula, is based on the Depth Domain Migration velocity model, is utilized three
Secondary Based on Interpolating Spline obtains the Migration velocity model of vertical time-domain, and then obtains the smoothed offset speed of vertical time-domain
Model.
Wherein the vertical time-domain coordinate conversion relation formula of Depth Domain-is：
In formula (1), x, y and z indicate the x, y and z directionss of cartesian coordinate system respectively；τ (x, y, z) is by Depth Domain coordinate
The vertical time-domain coordinate variable of variable z conversion gained, it is one-to-one with z；It is by the Depth Domain
Migration velocity model vdThe Depth Domain smoothed offset rate pattern of (x, y, z) smooth gained；Depth Domain can be sat using formula (1)
Mark is transformed to vertical time-domain；
It is described that vertical time-domain is obtained using cubic spline interpolation algorithm based on the Depth Domain Migration velocity model
Migration velocity model, and then the smoothed offset rate pattern of vertical time-domain is obtained, it specifically includes：
(a) according to the one-to-one relationship between Depth Domain coordinate z and vertical time-domain coordinate τ, the Depth Domain is smooth
Migration velocity modelThe Migration velocity model v of as original vertical time-domainτ1(x, y, τ), at this time vertical time
The coordinate interval of domain coordinate τ is non-uniform；
(b) the maximum coordinates interval for the vertical time-domain being subject to described in (a), using cubic spline interpolation algorithm to institute
The vertical time-domain coordinate stated is into row interpolation resampling, when obtaining uniform vertical time-domain coordinate, while can be obtained vertical
Between domain Migration velocity model vτ(x,y,τ)；
Similarly, by Depth Domain smoothed offset rate patternThe smoothed offset speed of vertical time-domain can be obtained
Spend model
Step 103：Isotropic medium one-order velocity-stress ACOUSTIC WAVE EQUATION is transformed into vertical time-domain by Depth Domain.
Using the vertical time-domain coordinate conversion relation formula of Depth Domain-, by isotropic medium one-order velocity-stress sound wave side of Depth Domain
Journey is transformed to isotropic medium one-order velocity-stress ACOUSTIC WAVE EQUATION of vertical time-domain.Wherein the Depth Domain is respectively to same
Property medium one-order velocity-stress ACOUSTIC WAVE EQUATION is specially：
In formula (2), ux, uyAnd uzRespectively Depth Domain isotropic medium sound wave particle velocity wave field is in x, y and z tri-
The component in direction, p are Depth Domain isotropic medium sound wave direct stress wave field；
It is described by isotropic medium one-order velocity-stress ACOUSTIC WAVE EQUATION of Depth Domain be transformed to vertical time-domain it is each to
Same sex medium one-order velocity-stress ACOUSTIC WAVE EQUATION, specially：
(a) Depth Domain coordinate system (x, y, z) is transformed to vertical time-domain coordinate systemThe seat of two coordinate systems
Mark variable meets following relationship：
(b) created symbol fxIt indicates：
Define vertical time-domain operator α=τxWith β=τy, according to chain rule, the partial derivative of cartesian coordinate system with it is vertical
The partial derivative of time-domain coordinate system meets following relationship：
So unit base vector of cartesian coordinate systemWith the unit base vector of bent coordinate systemBetween
Relationship is：
(c) the arbitrary vector under cartesian coordinate systemIt is in vertical time-domain corresponding form
For：
Wherein, Vx, VyAnd VzThe component in tri- directions x, y and z is tied up in cartesian coordinate for vector V；Vξ,And VτFor arrow
V is measured in vertical time-domain coordinate system ξ,With the component in tri- directions τ；
With divergence theorem, by the component of its vertical time-domain of the gradient of the vector V under cartesian coordinate system and divergence
And unit base vector is indicated, specially：
(d) by formula (8), (9) substitute into formula (2), you can obtain the vertical time domain one-order velocity-stress of isotropic medium
ACOUSTIC WAVE EQUATION, specially：
In formula (10), Uξ,UτRespectively vertical time-domain isotropic medium sound wave particle velocity wave field is when vertical
Between domain coordinate system ξ,The component in tri- directions τ, P are vertical time-domain isotropic medium sound wave direct stress wave field.
Step 104：Obtain speed-stress optimization operator per each corresponding moment of big gun.According to the observation system
Shot position coordinate of the parameter acquiring per big gun；Source wavelet is respectively set at each shot position coordinate；According to described vertical
Time-domain Migration velocity model, vertical time-domain smoothed offset rate pattern and offset parameter, utilize the preferential difference of staggered-mesh
Numerical algorithm solves the vertical time-domain single order stress-speed ACOUSTIC WAVE EQUATION, realizes the positive continuation per big gun source wavefield,
Obtain the source wavefield at each moment per big gun；Preserve the speed-stress optimization operator at each moment per big gun.It is wherein described
Arbitrary t moment, speed-stress optimization operator AtIt is defined as follows：
In formula (11),It is just answered for the vertical time-domain sound wave on i-th of mesh point outside t moment model area
Reeb；Wherein i=1,2 ..., N indicate that half exponent number of space difference used in wave field extrapolation, N indicate maximum half exponent number；H is indicated
Space lattice spacing；aiFor first group of Optimizing operator coefficient.
Step 105：Obtain the source wavefield at each moment of the big gun.Read the speed at each moment of every big gun of storage
Degree-stress optimization operator, the Optimizing operator using reading and specific difference scheme carry out source wavefield reconstruction, obtain every
The source wavefield of the reconstruction at each moment of big gun.
By taking direct stress wave field P as an example, the space in vertical time-domain isotropic medium single order stress-speed ACOUSTIC WAVE EQUATION
The staggering mesh finite-difference form of partial derivative is：By taking P is in x=0 (i.e. left margin) as an example, i-th point on the right side of x=0 of P
First-order partial derivative It can be indicated with finite difference, specially：
Wherein, bi,jFor second group of Optimizing operator coefficient, the as described specific difference scheme of formula (12), by formula (12)
With formula (11) comparison it is found that formula (12) contains the Optimizing operator At, wave field can be found out on a left side by formula (12) formula
The partial derivative of N-1 points in boundary；
Similarly, the partial derivative of N-1 points in model its coboundary can be obtained；
For the point inside model, its partial derivative is calculated using conventional staggering mesh finite-difference, to realize focus
Wave-field reconstruction obtains the source wavefield of the reconstruction at each moment per big gun.
Step 106：Obtain the geophone station wave field at each moment of the big gun.According to the observation system parameter acquiring per big gun
Geophone station position coordinates；According to the vertical time-domain Migration velocity model, vertical time-domain smoothed offset rate pattern and
Offset parameter solves the vertical time-domain single order stress-speed ACOUSTIC WAVE EQUATION, the sight to every big gun using numerical algorithm
Earthquake prediction record carries out inverse time continuation；Obtain the geophone station wave field at each moment per big gun.
Step 107：It obtains per the corresponding single-shot migrated section of big gun.Mutually in the same time, it is being based on migration imaging condition, to described
Every big gun each moment reconstruction source wavefield and geophone station wave field be imaged, obtain per big gun single-shot offset cuts open
Face；And then the single-shot migrated section of all big guns is obtained, the single-shot migrated section of all big guns is overlapped, vertical time-domain is obtained
Final migrated section.
The image-forming condition is：
In formula in (13), PS(x, y, τ, t) and PR(x, y, τ, t) indicates the arbitrary spatial position at each moment respectively
Vertical time-domain source wavefield and geophone station wave field, Is(x, y, τ) indicates that the single-shot migrated section of s big guns, T indicate dominant record
Time；
It is described to be overlapped the single-shot migrated section of all big guns, the final migrated section of vertical time-domain is obtained, is had
Body is：
In formula (14), I (x, y, τ) indicates the final migrated section of vertical time-domain.
Step 108：Obtain final migrated section.It, will be described using the vertical time-domain coordinate inverse transformation relationship of Depth Domain-
The final migration imaging section of vertical time-domain is converted into Depth Domain by vertical time-domain, obtains the final offset of Depth Domain
Section.
Wherein the vertical time-domain coordinate inverse transformation relationship of Depth Domain-is：
In formula (15), Z (x, y, τ) is by the Depth Domain coordinate variable of vertical time-domain coordinate τ inverse transformation Depth Domains；At this time
Z (x, y, τ) samples for unequal interval, and equal interval sampling is obtained into row interpolation to Z (x, y, τ) using cubic spline interpolation
Depth Domain coordinate variable
Using formula (15), by the final migration imaging section I (x, y, τ) of the vertical time-domain by vertical time-domain
It is converted into Depth Domain, obtains the final migrated section of Depth Domain
Specific implementation mode one：
Fig. 2 is the hollow Migration velocity model figure of two dimension that present embodiment provides, wherein Fig. 2 (a) is depth domain model,
Fig. 2 (b) is vertical time domain model.31 hypocenter of the explosions are set on this model, and source wavelet is set as Ricker wavelet, dominant frequency
It it is 30 hertz, starting focal point is located at (0m, 0m), and 200m is divided between big gun, and observation system is received using split shooting both sides, single
Side maximum offset is 3000m, and road spacing is 10m, and time step 1ms, record time span is 2s.
Fig. 3 is the single-shot earthquake record of the hollow model of two dimension shown in Fig. 2, wherein Fig. 3 (a) is the list obtained by conventional method
Big gun earthquake record, Fig. 3 (b) are the single-shot seismograms obtained by the method for the present invention.From figure 3, it can be seen that the method for the present invention can
To obtain the earthquake record with conventional High Resolution Finite Difference method same precision, it was demonstrated that the validity of 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 obtained by conventional method
Depth Domain migrated section, Fig. 4 (b) is the vertical time-domain migrated section obtained by the method for the present invention, and Fig. 4 (c) is present invention side
Depth Domain migrated section obtained by method.Figure 4, it is seen that the migrated section obtained by the method for the present invention and conventional method gained
Migrated section precision is identical, it was demonstrated that the method for the present invention has higher precision.In conjunction with earthquake record and migrated section comparison diagram
It is found that the method for the present invention can get high-precision migration imaging section while amount of storage is greatly reduced, there is very high meter
Calculate efficiency.
The present invention includes：Obtain seismic data, observation system parameter, the Migration velocity model of pending migration imaging
And offset parameter；Migration velocity model is transformed into vertical time-domain by Depth Domain；By isotropic medium one-order velocity-stress
ACOUSTIC WAVE EQUATION transforms to vertical time-domain by Depth Domain；Obtain speed-stress optimization operator per each corresponding moment of big gun；
Source wavefield reconstruction is carried out using the speed-stress optimization operator at each corresponding moment of the big gun of gained, obtains the big gun
The source wavefield at each moment；Obtain the geophone station wave field at each moment per big gun；Obtain the corresponding single-shot offset per big gun
Section；The single-shot migrated section of all big guns is overlapped, final migrated section is obtained.The method of the present invention, Ke Yiyou
Effect reduces the calculation amount and amount of storage of conventional three-dimensional reverse-time migration method, improves computational efficiency.
Each embodiment is described by the way of progressive in this specification, the highlights of each of the examples are with other
The difference of embodiment, just to refer each other for identical similar portion between each embodiment.For system disclosed in embodiment
For, since it is corresponded to the methods disclosed in the examples, so description is fairly simple, related place is said referring to method part
It is bright.
Principle and implementation of the present invention are described for specific case used herein, and above example is said
The bright method and its core concept for being merely used to help understand the present invention；Meanwhile for those of ordinary skill in the art, foundation
The thought of the present invention, there will be changes in the specific implementation manner and application range.In conclusion the content of the present specification is not
It is interpreted as limitation of the present invention.
Claims (7)
1. a kind of seismic data offset imaging method, which is characterized in that the method includes：
Read earthquake work area observation system parameter, rate pattern and more big gun seismic datas of observation, determine into line displacement at
More big guns observation earthquake record, observation system parameter, Depth Domain Migration velocity model and the offset parameter of picture；
Using the vertical time-domain coordinate conversion relation formula of Depth Domain-, Depth Domain coordinate is transformed to vertical time-domain, based on described
Depth Domain Migration velocity model obtains the Migration velocity model of vertical time-domain, and then obtain using cubic spline interpolation algorithm
The smoothed offset rate pattern of vertical time-domain；
Using the vertical time-domain coordinate conversion relation formula of Depth Domain-, by isotropic medium one-order velocity-stress sound of Depth Domain
Wave equation is transformed to isotropic medium one-order velocity-stress ACOUSTIC WAVE EQUATION of vertical time-domain；
Shot position coordinate according to the observation system parameter acquiring per big gun；Shake is respectively set at each shot position coordinate
Source wavelet；According to the vertical time-domain Migration velocity model, vertical time-domain smoothed offset rate pattern and offset parameter, profit
The vertical time-domain single order stress-speed ACOUSTIC WAVE EQUATION is solved with the preferential difference numerical algorithm of staggered-mesh, is realized per big gun
The positive continuation of source wavefield obtains the source wavefield at each moment per big gun；The speed-for preserving each moment per big gun is answered
Power Optimizing operator；
Speed-stress optimization the operator for reading each moment of every big gun of storage utilizes the Optimizing operator of reading and specific
Difference scheme, carry out source wavefield reconstruction, obtain per big gun each moment reconstruction source wavefield；
Geophone station position coordinates according to the observation system parameter acquiring per big gun；According to the vertical time-domain migration velocity mould
Type, vertical time-domain smoothed offset rate pattern and offset parameter solve the vertical time-domain single order using numerical algorithm
Stress-speed ACOUSTIC WAVE EQUATION carries out inverse time continuation to the observation earthquake record of every big gun；Obtain each moment per big gun
Geophone station wave field；
Mutually in the same time, it is being based on migration imaging condition, the source wavefield to the reconstruction at each moment of every big gun and inspection
Wave point wave field is imaged, and the single-shot migrated section per big gun is obtained；And then the single-shot migrated section of all big guns is obtained, by all big guns
Single-shot migrated section be overlapped, obtain the final migrated section of vertical time-domain；
Using the vertical time-domain coordinate inverse transformation relationship of Depth Domain-, by the final migration imaging section of the vertical time-domain
Depth Domain is converted by vertical time-domain, obtains the final migrated section of Depth Domain.
2. a kind of seismic data offset imaging method according to claim 1, which is characterized in that described to be hung down using Depth Domain-
To time-domain coordinate conversion relation formula, Depth Domain coordinate is transformed to vertical time-domain, is based on the Depth Domain migration velocity mould
Type obtains the Migration velocity model of vertical time-domain, and then obtain the smooth of vertical time-domain using cubic spline interpolation algorithm
Migration velocity model, specially：
The vertical time-domain coordinate conversion relation formula of Depth Domain-is：
In formula (1), x, y and z indicate the x, y and z directionss of cartesian coordinate system respectively；τ (x, y, z) is by Depth Domain coordinate variable
The vertical time-domain coordinate variable of z conversion gained, it is one-to-one with z；It is to be deviated by the Depth Domain
Rate pattern vdThe Depth Domain smoothed offset rate pattern of (x, y, z) smooth gained；Depth Domain coordinate can be become using formula (1)
It is changed to vertical time-domain；
It is described that the offset of vertical time-domain is obtained using cubic spline interpolation algorithm based on the Depth Domain Migration velocity model
Rate pattern, and then the smoothed offset rate pattern of vertical time-domain is obtained, it specifically includes：
(a) according to the one-to-one relationship between Depth Domain coordinate z and vertical time-domain coordinate τ, the Depth Domain smoothed offset
Rate patternThe Migration velocity model v of as original vertical time-domainτ1(x, y, τ), at this time vertical time-domain sit
The coordinate interval for marking τ is non-uniform；
(b) the maximum coordinates interval for the vertical time-domain being subject to described in (a), using cubic spline interpolation algorithm to described
Vertical time-domain coordinate obtains uniform vertical time-domain coordinate, while can be obtained vertical time-domain into row interpolation resampling
Migration velocity model vτ(x,y,τ)；
Similarly, by Depth Domain smoothed offset rate patternThe smoothed offset speed mould of vertical time-domain can be obtained
Type
3. according to a kind of seismic data offset imaging method described in claim 1, it is characterised in that：It is described to be hung down using Depth Domain-
To time-domain coordinate conversion relation formula, isotropic medium one-order velocity-stress ACOUSTIC WAVE EQUATION of Depth Domain is transformed to vertical
Isotropic medium one-order velocity-stress ACOUSTIC WAVE EQUATION of time-domain, specially：
Depth Domain isotropic medium one-order velocity-stress ACOUSTIC WAVE EQUATION is specially：
In formula (2), ux, uyAnd uzRespectively Depth Domain isotropic medium sound wave particle velocity wave field is in tri- directions x, y and z
Component, p be Depth Domain isotropic medium sound wave direct stress wave field；
The isotropism that isotropic medium one-order velocity-stress ACOUSTIC WAVE EQUATION of Depth Domain is transformed to vertical time-domain
Medium one-order velocity-stress ACOUSTIC WAVE EQUATION, specially：
(a) Depth Domain coordinate system (x, y, z) is transformed to vertical time-domain coordinate systemThe coordinate variable of two coordinate systems
Meet following relationship：
(b) created symbol fxIt indicates：
Define vertical time-domain operator α=τxWith β=τy, according to chain rule, the partial derivative of cartesian coordinate system and vertical time
The partial derivative of domain coordinate system meets following relationship：
So unit base vector of cartesian coordinate systemWith the unit base vector of bent coordinate systemBetween relationship
For：
(c) the arbitrary vector under cartesian coordinate systemIt is in vertical time-domain corresponding form：
Wherein, Vx, VyAnd VzThe component in tri- directions x, y and z is tied up in cartesian coordinate for vector V；Vξ,And VτFor vector V
In vertical time-domain coordinate system ξ,With the component in tri- directions τ；
With divergence theorem, by the component and list of its vertical time-domain of the gradient of the vector V under cartesian coordinate system and divergence
Position base vector is indicated, specially：
(d) by formula (8), (9) substitute into formula (2), you can obtain the vertical time domain one-order velocity-stress sound wave of isotropic medium
Equation, specially：
In formula (10), Uξ,UτRespectively vertical time-domain isotropic medium sound wave particle velocity wave field is in vertical time-domain
Coordinate system ξ,The component in tri- directions τ, P are vertical time-domain isotropic medium sound wave direct stress wave field.
4. according to a kind of seismic data offset imaging method described in claim 1, it is characterised in that：It is described to preserve each of every big gun
The speed at moment-stress optimization operator, specially：
Arbitrary t moment, speed-stress optimization operator AtIt is defined as follows：
In formula (11),For the vertical time-domain sound wave direct stress wave on i-th of mesh point outside t moment model area
；Wherein i=1,2 ..., N indicate that half exponent number of space difference used in wave field extrapolation, N indicate maximum half exponent number；H representation spaces
Grid spacing；aiFor first group of Optimizing operator coefficient.
5. according to a kind of seismic data offset imaging method described in claim 1, it is characterised in that：It is described to utilize the optimization read
Operator and specific difference scheme carry out source wavefield reconstruction, obtain the source wavefield of the reconstruction at each moment per big gun,
Specially；
By taking direct stress wave field P as an example, the space local derviation in vertical time-domain isotropic medium single order stress-speed ACOUSTIC WAVE EQUATION
Several staggering mesh finite-difference forms are：By taking P is in x=0 (i.e. left margin) as an example, i-th point of single orders of the P on the right side of x=0
Partial derivative It can be indicated with finite difference, specially：
Wherein, bi,jFor second group of Optimizing operator coefficient, the as described specific difference scheme of formula (12), by formula (12) and formula
Sub (11) comparison is it is found that formula (12) contains the Optimizing operator At, wave field can be found out in left margin by formula (12) formula
The partial derivative of interior N-1 points；
Similarly, the partial derivative of N-1 points in model its coboundary can be obtained；
For the point inside model, its partial derivative is calculated using conventional staggering mesh finite-difference, to realize source wavefield
It rebuilds, obtains the source wavefield of the reconstruction at each moment per big gun.
6. according to a kind of seismic data offset imaging method described in claim 1, it is characterised in that：It is described mutually in the same time, base
In migration imaging condition, the source wavefield and geophone station wave field of the reconstruction at each moment of every big gun are imaged,
Obtain the single-shot migrated section per big gun；And then obtain the single-shot migrated sections of all big guns, by the single-shot migrated section of all big guns into
Row superposition, obtains the final migrated section of vertical time-domain, specially：
The image-forming condition is：
In formula in (13), PS(x, y, τ, t) and PR(x, y, τ, t) indicates the vertical of the arbitrary spatial position at each moment respectively
Time-domain source wavefield and geophone station wave field, Is(x, y, τ) indicates the single-shot migrated section of s big guns, when T indicates dominant record
Between；
It is described to be overlapped the single-shot migrated section of all big guns, the final migrated section of vertical time-domain is obtained, specially：
In formula (14), I (x, y, τ) indicates the final migrated section of vertical time-domain.
7. according to a kind of seismic data offset imaging method described in claim 1, it is characterised in that：It is described to be hung down using Depth Domain-
To time-domain coordinate inverse transformation relationship, the final migration imaging section of the vertical time-domain is converted by vertical time-domain
Depth Domain obtains the final migrated section of Depth Domain, specially：
The vertical time-domain coordinate inverse transformation relationship of Depth Domain-is：
In formula (15), Z (x, y, τ) is by the Depth Domain coordinate variable of vertical time-domain coordinate τ inverse transformation Depth Domains；At this time Z (x,
Y, τ) it is that unequal interval samples, the depth of equal interval sampling is obtained into row interpolation to Z (x, y, τ) using cubic spline interpolation
Domain coordinate variable
Using formula (15), the final migration imaging section I (x, y, τ) of the vertical time-domain is converted by vertical time-domain
To Depth Domain, the final migrated section of Depth Domain is obtained
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810104681.4A CN108363097B (en) | 2018-02-02 | 2018-02-02 | A kind of seismic data offset imaging method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810104681.4A CN108363097B (en) | 2018-02-02 | 2018-02-02 | A kind of seismic data offset imaging method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108363097A true CN108363097A (en) | 2018-08-03 |
| CN108363097B CN108363097B (en) | 2019-10-15 |
Family
ID=63004198
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810104681.4A Active CN108363097B (en) | 2018-02-02 | 2018-02-02 | A kind of seismic data offset imaging method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108363097B (en) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111158049A (en) * | 2019-12-27 | 2020-05-15 | 同济大学 | Seismic reverse time migration imaging method based on scattering integration method |
| CN111337973A (en) * | 2018-12-19 | 2020-06-26 | 中国石油天然气集团有限公司 | Seismic data reconstruction method and system |
| CN112946736A (en) * | 2019-11-26 | 2021-06-11 | 中国石油天然气集团有限公司 | Three-dimensional observation system reconstruction method and system |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| FR2666905A1 (en) * | 1990-09-18 | 1992-03-20 | Total Petroles | Improved method of underground seismic prospecting |
| US20040196738A1 (en) * | 2003-04-07 | 2004-10-07 | Hillel Tal-Ezer | Wave migration by a krylov space expansion of the square root exponent operator, for use in seismic imaging |
| CN105652321A (en) * | 2015-12-30 | 2016-06-08 | 中国石油大学(华东) | Visco-acoustic anisotropic least square inverse time migration imaging method |
| CN106033124A (en) * | 2016-06-29 | 2016-10-19 | 中国石油化工股份有限公司 | Multi-seismic resource sticky sound least square reverse time migration method based on stochastic optimization |
| CN106970416A (en) * | 2017-03-17 | 2017-07-21 | 中国地质科学院地球物理地球化学勘查研究所 | Elastic wave least square reverse-time migration system and method based on wave field separation |
-
2018
- 2018-02-02 CN CN201810104681.4A patent/CN108363097B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| FR2666905A1 (en) * | 1990-09-18 | 1992-03-20 | Total Petroles | Improved method of underground seismic prospecting |
| US20040196738A1 (en) * | 2003-04-07 | 2004-10-07 | Hillel Tal-Ezer | Wave migration by a krylov space expansion of the square root exponent operator, for use in seismic imaging |
| CN105652321A (en) * | 2015-12-30 | 2016-06-08 | 中国石油大学(华东) | Visco-acoustic anisotropic least square inverse time migration imaging method |
| CN106033124A (en) * | 2016-06-29 | 2016-10-19 | 中国石油化工股份有限公司 | Multi-seismic resource sticky sound least square reverse time migration method based on stochastic optimization |
| CN106970416A (en) * | 2017-03-17 | 2017-07-21 | 中国地质科学院地球物理地球化学勘查研究所 | Elastic wave least square reverse-time migration system and method based on wave field separation |
Non-Patent Citations (1)
| Title |
|---|
| 杨凯 等: "基于非结构化网格的最小二乘逆时偏移", 《地球物理学报》 * |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111337973A (en) * | 2018-12-19 | 2020-06-26 | 中国石油天然气集团有限公司 | Seismic data reconstruction method and system |
| CN111337973B (en) * | 2018-12-19 | 2022-08-05 | 中国石油天然气集团有限公司 | Seismic data reconstruction method and system |
| CN112946736A (en) * | 2019-11-26 | 2021-06-11 | 中国石油天然气集团有限公司 | Three-dimensional observation system reconstruction method and system |
| CN111158049A (en) * | 2019-12-27 | 2020-05-15 | 同济大学 | Seismic reverse time migration imaging method based on scattering integration method |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108363097B (en) | 2019-10-15 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107783190B (en) | A kind of least square reverse-time migration gradient updating method | |
| CN106526674B (en) | Three-dimensional full waveform inversion energy weighting gradient preprocessing method | |
| CN108108331B (en) | A kind of finite difference formulations method based on quasi- spatial domain equations for elastic waves | |
| CN108241173B (en) | A kind of seismic data offset imaging method and system | |
| US8103453B2 (en) | Method of seismic data interpolation by projection on convex sets | |
| CN108051855B (en) | A kind of finite difference formulations method based on quasi- spatial domain ACOUSTIC WAVE EQUATION | |
| CN108363097B (en) | A kind of seismic data offset imaging method | |
| CN104459770B (en) | A kind of method for regularizing high-dimensional seismic data | |
| CN105607122B (en) | A kind of earthquake texture blending and Enhancement Method based on full variation geological data decomposition model | |
| WO2017195093A1 (en) | Simultaneous source acquisition and separation on general related sampling grids | |
| CN110954945B (en) | Full waveform inversion method based on dynamic random seismic source coding | |
| CN110579795B (en) | Joint velocity inversion method based on passive source seismic waveform and reverse-time imaging thereof | |
| CN106932819A (en) | Pre-stack seismic parameter inversion method based on anisotropy Markov random field | |
| WO2011103297A2 (en) | Estimating internal multiples in seismic data | |
| CN109541681A (en) | A kind of waveform inversion method of streamer seismic data and a small amount of OBS data aggregate | |
| CN102830431B (en) | Self-adaption interpolating method for real ground-surface ray tracking | |
| CN109633749A (en) | Non-linear Fresnel zone seismic traveltime tomography method based on scattering integral method | |
| CN107356972A (en) | A kind of imaging method of anisotropic medium | |
| CN113740915B (en) | Method for jointly inverting crust structure parameters by spherical coordinate system gravity and receiving function | |
| Hole et al. | Interface inversion using broadside seismic refraction data and three‐dimensional travel time calculations | |
| Xiang‐Bo et al. | Anisotropic Radon transform and its application to demultiple | |
| CN106842300A (en) | A kind of high efficiency multi-component seismic data true amplitude migration imaging method | |
| CN107340537A (en) | A kind of method of P-SV converted waves prestack reverse-time depth migration | |
| Qu et al. | 3-D Least-Squares Reverse Time Migration in Curvilinear-τ Domain | |
| CN109557581A (en) | Reconstruction of seismic data method and system based on Fourier transformation |
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 |