CN108363097B - A kind of seismic data offset imaging method - Google Patents
A kind of seismic data offset imaging method Download PDFInfo
- Publication number
- CN108363097B CN108363097B CN201810104681.4A CN201810104681A CN108363097B CN 108363097 B CN108363097 B CN 108363097B CN 201810104681 A CN201810104681 A CN 201810104681A CN 108363097 B CN108363097 B CN 108363097B
- Authority
- CN
- China
- Prior art keywords
- domain
- vertical time
- depth
- coordinate
- time
- 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.)
- Active
Links
- 238000003384 imaging method Methods 0.000 title claims abstract description 30
- 238000013508 migration Methods 0.000 claims abstract description 67
- 230000005012 migration Effects 0.000 claims abstract description 67
- 238000000034 method Methods 0.000 claims abstract description 37
- 238000005457 optimization Methods 0.000 claims abstract description 17
- 238000006243 chemical reaction Methods 0.000 claims description 14
- 230000009466 transformation Effects 0.000 claims description 11
- 239000002245 particle Substances 0.000 claims description 6
- 238000005070 sampling Methods 0.000 claims description 5
- 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
- 238000007689 inspection Methods 0.000 claims 2
- 238000004364 calculation method Methods 0.000 abstract description 4
- 238000005516 engineering process Methods 0.000 description 7
- 238000007796 conventional method Methods 0.000 description 5
- 230000008901 benefit Effects 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 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
Images
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 comprises: obtaining the seismic data of pending migration imaging, observation system parameter, Migration velocity model and offset parameter；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 at each corresponding moment of every big gun；Source wavefield reconstruction is carried out using the speed-stress optimization operator at each corresponding moment of the resulting big gun, obtains the source wavefield at each moment of the big gun；Obtain the geophone station wave field at each moment of every big gun；Obtain the corresponding single-shot migrated section of every big gun；The single-shot migrated section of all big guns is overlapped, final migrated section is obtained.The calculation amount and amount of storage of conventional three-dimensional reverse-time migration method can be effectively reduced in the method for the present invention, improves and calculates efficiency.
Description
Technical field
The invention belongs to seismic exploration technique fields, more particularly to a kind of seismic data offset imaging method.
Background technique
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, at present common several method technology in the industry are as follows: Checkpoint technology, RANDOM BOUNDARY wave-field reconstruction skill
Art, efficiency frontier memory technology.These types of technology has respective limitation: Checkpoint technology calculation amount is 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 amount of storage present in current reverse-time migration method and computational efficiency problem, 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 ".
Summary of the invention
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, which comprises
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；
The shot position coordinate of every big gun is obtained according to the observation system parameter；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
The positive continuation of every big gun source wavefield, obtains the source wavefield at each moment of every big gun；Save the speed at each moment of every big gun
Degree-stress optimization operator；
The geophone station position coordinates of every big gun are obtained according to the observation system parameter；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 of every 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 of every 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 of every 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, specifically:
The vertical time-domain coordinate conversion relation formula of Depth Domain-are as follows:
In formula (1), x, y and z respectively indicate the x, y and z directionss of cartesian coordinate system；τ (x, y, z) is by Depth Domain coordinate
Variable z converts resulting vertical time-domain coordinate variable, it is one-to-one with z；
It is by the depth
Domain Migration velocity model vd(x, y, z) smooth resulting Depth Domain smoothed offset rate pattern；It can be by Depth Domain using formula (1)
Coordinate 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 model
The 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 of vertical time-domain described in (a) is subject to, using cubic spline interpolation algorithm to institute
The vertical time-domain coordinate stated carries out 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 pattern
The 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
Are as follows:
Depth Domain isotropic medium one-order velocity-stress ACOUSTIC WAVE EQUATION specifically:
In formula (2), ux, uyAnd uzRespectively Depth Domain isotropic medium sound wave particle velocity wave field is at 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, specifically:
(a) Depth Domain coordinate system (x, y, z) is transformed to vertical time-domain coordinate system
The 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 system
With the unit base vector of bent coordinate system
Between
Relationship are as follows:
(c) any vector under cartesian coordinate system
It is in vertical time-domain corresponding form
Are as follows:
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, specifically:
(d) by formula (8), (9) substitute into formula (2), the vertical time domain one-order velocity-stress of isotropic medium can be obtained
ACOUSTIC WAVE EQUATION, specifically:
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 at each moment for saving every big gun, specifically:
Any 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 difference half order in space used in wave field extrapolation, N indicate maximum half order；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 of every big gun, specially；
Space by taking direct stress wave field P as an example, in vertical time-domain isotropic medium single order stress-speed ACOUSTIC WAVE EQUATION
The staggering mesh finite-difference form of partial derivative are as follows: by taking P is at 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, specifically:
Wherein, bi,jFor second group of Optimizing operator coefficient, formula (12) is the specific difference scheme, 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 point in boundary；
Similarly, the partial derivative of N-1 point in its coboundary of model 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 of every 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 the single-shot migrated section of every big gun；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, the final migrated section of vertical time-domain is obtained, specifically:
The image-forming condition are as follows:
In formula in (13), PS(x, y, τ, t) and PR(x, y, τ, t) respectively indicates any spatial position at each moment
Vertical time-domain source wavefield and geophone station wave field, Is(x, y, τ) indicates that the single-shot migrated section of s big gun, 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 are as follows:
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, specifically:
The vertical time-domain coordinate inverse transformation relationship of Depth Domain-are as follows:
In formula (15), Z (x, y, τ) is by the Depth Domain coordinate variable of vertical time-domain coordinate τ inverse transformation Depth Domain；At this time
Z (x, y, τ) is unequal interval sampling, carries out interpolation to Z (x, y, τ) using cubic spline interpolation, obtains equal interval sampling
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 technical 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 benefit
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, can get the migration imaging section of high quality；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.
Detailed description of the invention
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 invention
Example, for those of ordinary skill in the art, without any 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
It (b) 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 resulting list of conventional method
Big gun earthquake record, Fig. 3 (b) are the resulting single-shot seismograms 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 resulting vertical time-domain migrated section of the method for the present invention, and Fig. 4 (c) is present invention side
The resulting Depth Domain migrated section of method.
Specific embodiment
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 description, it is clear that described embodiment is only a part of the embodiments 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 with reference to the accompanying drawing and specific real to keep above-mentioned purpose of the invention, features and advantages 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, which comprises
Step 101: reading earthquake work area observation system parameter, rate pattern and more big gun seismic datas of observation, 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.
The wherein vertical time-domain coordinate conversion relation formula of Depth Domain-are as follows:
In formula (1), x, y and z respectively indicate the x, y and z directionss of cartesian coordinate system；τ (x, y, z) is by Depth Domain coordinate
Variable z converts resulting vertical time-domain coordinate variable, it is one-to-one with z；
It is by the depth
Domain Migration velocity model vd(x, y, z) smooth resulting Depth Domain smoothed offset rate pattern；It can be by Depth Domain using formula (1)
Coordinate 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 model
The 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 of vertical time-domain described in (a) is subject to, using cubic spline interpolation algorithm to institute
The vertical time-domain coordinate stated carries out 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 pattern
The 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 specifically:
In formula (2), ux, uyAnd uzRespectively Depth Domain isotropic medium sound wave particle velocity wave field is at 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, specifically:
(a) Depth Domain coordinate system (x, y, z) is transformed to vertical time-domain coordinate system
The 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 system
With the unit base vector of bent coordinate system
Between
Relationship are as follows:
(c) any vector under cartesian coordinate system
It is in vertical time-domain corresponding form
Are as follows:
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 of its vertical time-domain of the gradient of the vector V under cartesian coordinate system and divergence
And unit base vector is indicated, specifically:
(d) by formula (8), (9) substitute into formula (2), the vertical time domain one-order velocity-stress of isotropic medium can be obtained
ACOUSTIC WAVE EQUATION, specifically:
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: obtaining the speed-stress optimization operator at each corresponding moment of every big gun.According to the observation system
Parameter obtains the shot position coordinate of every 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 of every big gun source wavefield,
Obtain the source wavefield at each moment of every big gun；Save the speed-stress optimization operator at each moment of every big gun.It is wherein described
Any 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 difference half order in space used in wave field extrapolation, N indicate maximum half order；H is indicated
Space lattice spacing；aiFor first group of Optimizing operator coefficient.
Step 105: obtaining 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, Optimizing operator and specific difference scheme using reading carry out source wavefield reconstruction, obtain every
The source wavefield of the reconstruction at each moment of big gun.
Space by taking direct stress wave field P as an example, in vertical time-domain isotropic medium single order stress-speed ACOUSTIC WAVE EQUATION
The staggering mesh finite-difference form of partial derivative are as follows: by taking P is at 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, specifically:
Wherein, bi,jFor second group of Optimizing operator coefficient, formula (12) is the specific difference scheme, 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 point in boundary；
Similarly, the partial derivative of N-1 point in its coboundary of model 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 of every big gun.
Step 106: obtaining the geophone station wave field at each moment of the big gun.Every big gun is obtained according to the observation system parameter
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 of every big gun.
Step 107: obtaining the corresponding single-shot migrated section of every 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 every 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 are as follows:
In formula in (13), PS(x, y, τ, t) and PR(x, y, τ, t) respectively indicates any spatial position at each moment
Vertical time-domain source wavefield and geophone station wave field, Is(x, y, τ) indicates that the single-shot migrated section of s big gun, 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 are as follows:
In formula (14), I (x, y, τ) indicates the final migrated section of vertical time-domain.
Step 108: obtaining 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.
The wherein vertical time-domain coordinate inverse transformation relationship of Depth Domain-are as follows:
In formula (15), Z (x, y, τ) is by the Depth Domain coordinate variable of vertical time-domain coordinate τ inverse transformation Depth Domain；At this time
Z (x, y, τ) is unequal interval sampling, carries out interpolation to Z (x, y, τ) using cubic spline interpolation, obtains equal interval sampling
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 embodiment 1:
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 is 30 hertz, starting focal point is located at (0m, 0m), and 200m is divided between big gun, receives observation system 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 resulting list of conventional method
Big gun earthquake record, Fig. 3 (b) are the resulting single-shot seismograms of 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 resulting vertical time-domain migrated section of the method for the present invention, and Fig. 4 (c) is present invention side
The resulting Depth Domain migrated section of method.Figure 4, it is seen that obtained by the resulting migrated section of the method for the present invention and conventional method
Migrated section precision is identical, it was demonstrated that the method for the present invention precision with higher.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: the seismic data for obtaining pending migration imaging, observation system parameter, Migration velocity model
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 the speed-stress optimization operator at each corresponding moment of every big gun；
Source wavefield reconstruction is carried out using the speed-stress optimization operator at each corresponding moment of the resulting big gun, obtains the big gun
The source wavefield at each moment；Obtain the geophone station wave field at each moment of every big gun；Obtain the corresponding single-shot offset of every 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 in this specification is described in a progressive manner, the highlights of each of the examples are with other
The difference of embodiment, the same or similar parts in each embodiment may refer to each other.For system disclosed in embodiment
For, since it is corresponded to the methods disclosed in the examples, so being described relatively simple, related place is said referring to method part
It is bright.
Used herein a specific example illustrates the principle and implementation of the invention, and above embodiments are said
It is bright to be merely used to help understand method and its core concept of the invention；At the same time, for those skilled in the art, foundation
Thought of the 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 (1)
1. a kind of seismic data offset imaging method, which is characterized in that the described method includes:
Earthquake work area observation system parameter, rate pattern and more big gun seismic datas of observation are read, determination is deviated 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 coordinate, is based on
The Depth Domain Migration velocity model obtains vertical time-domain Migration velocity model using cubic spline interpolation algorithm, and then obtains
Obtain vertical time-domain smoothed offset rate pattern；
Using the vertical time-domain coordinate conversion relation formula of Depth Domain-, by Depth Domain isotropic medium one-order velocity-stress sound wave
Equation transform is vertical time-domain isotropic medium one-order velocity-stress ACOUSTIC WAVE EQUATION；
The shot position coordinate of every big gun is obtained according to the observation system parameter；Shake is respectively set at each shot position coordinate
Source wavelet；According to vertical time-domain Migration velocity model, vertical time-domain smoothed offset rate pattern and offset parameter, friendship is utilized
The wrong preferential difference numerical algorithm of grid solves vertical time-domain single order stress-speed ACOUSTIC WAVE EQUATION, realizes every big gun source wavefield
Positive continuation obtains the source wavefield at each moment of every big gun；Speed-the stress optimization for saving each moment of every big gun is calculated
Son；
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 the source wavefield of the reconstruction at each moment of every big gun；
The geophone station position coordinates of every big gun are obtained according to the observation system parameter；According to vertical time-domain Migration velocity model,
Vertical time-domain smoothed offset rate pattern and offset parameter solve vertical time-domain single order stress-speed using numerical algorithm
ACOUSTIC WAVE EQUATION carries out inverse time continuation to the observation earthquake record of every big gun；Obtain the geophone station wave field at each moment of every big gun；
Mutually in the same time, it is being based on migration imaging condition, source wavefield and inspection to the reconstruction at each moment of every big gun
Wave point wave field is imaged, and the single-shot migrated section of every 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 migrated section of vertical time-domain by it is vertical when
Between domain be converted into Depth Domain, obtain the final migrated section of Depth Domain；
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 vertical time-domain Migration velocity model using cubic spline interpolation algorithm, and then is hung down
To time-domain smoothed offset rate pattern, comprising:
The vertical time-domain coordinate conversion relation formula of Depth Domain-are as follows:
In formula (1), x, y and z respectively indicate the x, y and z directionss of cartesian coordinate system；τ (x, y, z) is by Depth Domain coordinate variable
Convert resulting vertical time-domain coordinate variable；
It is by the Depth Domain Migration velocity model vd(x,y,z)
Smooth resulting Depth Domain smoothed offset rate pattern；Depth Domain coordinate can be transformed to vertical time-domain using formula (1)；
It is described to be based on the Depth Domain Migration velocity model, using cubic spline interpolation algorithm, obtain vertical time-domain offset speed
Model is spent, and then obtains vertical time-domain smoothed offset rate pattern, is specifically included:
(a) according to the one-to-one relationship between Depth Domain coordinate and vertical time-domain coordinate, the Depth Domain smoothed offset speed
Spend model
As original vertical time-domain Migration velocity model vτ1(x, y, τ), at this time vertical time-domain coordinate
Coordinate interval is non-uniform；
(b) the maximum coordinates interval of vertical time-domain in step (a) is subject to, using cubic spline interpolation algorithm to the vertical time
Domain coordinate carries out interpolation resampling, obtains uniform vertical time-domain coordinate, while can be obtained vertical time-domain migration velocity
Model vτ(x,y,τ)；
Similarly, by Depth Domain smoothed offset rate pattern
Vertical time-domain smoothed offset rate pattern can be obtained
Using the vertical time-domain coordinate conversion relation formula of Depth Domain-, by Depth Domain isotropic medium one-order velocity-stress sound wave
Equation transform is vertical time-domain isotropic medium one-order velocity-stress ACOUSTIC WAVE EQUATION, specifically:
Depth Domain isotropic medium one-order velocity-stress ACOUSTIC WAVE EQUATION are as follows:
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；
It is described that Depth Domain isotropic medium one-order velocity-stress ACOUSTIC WAVE EQUATION is transformed to vertical time-domain isotropic medium
One-order velocity-stress ACOUSTIC WAVE EQUATION, specifically:
(a) Depth Domain coordinate system (x, y, z) is transformed to vertical time-domain coordinate system
The coordinate of two coordinate systems becomes
Amount 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 and vertical time
The partial derivative of domain coordinate system meets following relationship:
So unit base vector of cartesian coordinate system
With the unit base vector of bent coordinate system
Between relationship
Are as follows:
(c) any vector under cartesian coordinate systemIt is in vertical time-domain corresponding form are as follows:
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, specifically:
(d) by formula (8), (9) substitute into formula (2), the vertical time domain one-order velocity-stress sound wave of isotropic medium can be obtained
Equation, specifically:
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；
Speed-stress optimization the operator at each moment of every big gun is saved, specifically:
Any 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 difference half order in space used in wave field extrapolation, N indicate maximum half order；H representation space
Grid spacing；aiFor first group of Optimizing operator coefficient；
Optimizing operator and specific difference scheme using reading carry out source wavefield reconstruction, when obtaining each of every big gun
The source wavefield of the reconstruction at quarter, specially；
Space local derviation when direct stress wave field is P, in vertical time-domain isotropic medium single order stress-speed ACOUSTIC WAVE EQUATION
Several staggering mesh finite-difference form are as follows: P is in x=0, then i-th point of first-order partial derivative of the P on the right side of x=0
It can be indicated with finite difference, specifically:
Wherein, bi,jFor second group of Optimizing operator coefficient, formula (12) is specific difference scheme, by formula (12) and formula (11)
Comparison is it is found that formula (12) contains the Optimizing operator At, wave field N-1 point in left margin can be found out by formula (12) formula
Partial derivative；
Similarly, the partial derivative of N-1 point in its coboundary of model 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 of every big gun；
Mutually in the same time, it is being based on migration imaging condition, source wavefield and inspection to the reconstruction at each moment of every big gun
Wave point wave field is imaged, and the single-shot migrated section of every 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, specifically:
The migration imaging condition are as follows:
In formula in (13), PS(x, y, τ, t) and PR(x, y, τ, t) respectively indicates the vertical of any spatial position at each moment
Time-domain source wavefield and geophone station wave field, Is(x, y, τ) indicates the single-shot migrated section of s big gun, 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, specifically:
In formula (14), I (x, y, τ) indicates the final migrated section of vertical time-domain；
Using the vertical time-domain coordinate inverse transformation relationship of Depth Domain-, by the final migrated section of vertical time-domain by it is vertical when
Between domain be converted into Depth Domain, obtain the final migrated section of Depth Domain, specifically:
The vertical time-domain coordinate inverse transformation relationship of Depth Domain-are as follows:
In formula (15), Z (x, y, τ) is by the Depth Domain coordinate variable of vertical time-domain coordinate τ inverse transformation Depth Domain；At this time Z (x,
Y, τ) it is that unequal interval samples, interpolation is carried out to Z (x, y, τ) using cubic spline interpolation, obtains the depth of equal interval sampling
Domain coordinate variable
Using formula (15), the final migrated section I (x, y, τ) of vertical time-domain is converted into depth by vertical time-domain
Domain obtains the final migrated section of Depth Domain
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 CN108363097A (en) | 2018-08-03 |
| CN108363097B true 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) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| 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 |
| CN111158049B (en) * | 2019-12-27 | 2020-11-27 | 同济大学 | Seismic reverse time migration imaging method based on scattering integration method |
Family Cites Families (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| FR2666905B1 (en) * | 1990-09-18 | 1992-12-24 | Total Petroles | IMPROVED SEISMIC PROSPECTION PROCESS OF THE BASEMENT. |
| 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 |
| CN105652321B (en) * | 2015-12-30 | 2016-10-12 | 中国石油大学(华东) | A kind of viscous acoustic anisotropy least square reverse-time migration formation method |
| CN106033124B (en) * | 2016-06-29 | 2018-10-16 | 中国石油化工股份有限公司 | A kind of viscous sound least square reverse-time migration method of more focus based on stochastic optimization |
| CN106970416B (en) * | 2017-03-17 | 2018-12-04 | 中国地质科学院地球物理地球化学勘查研究所 | 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
Also Published As
| Publication number | Publication date |
|---|---|
| CN108363097A (en) | 2018-08-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Hennenfent et al. | Nonequispaced curvelet transform for seismic data reconstruction: A sparsity-promoting approach | |
| CN108241173B (en) | A kind of seismic data offset imaging method and system | |
| CN103091710B (en) | A kind of reverse-time migration formation method and device | |
| CN106526677B (en) | A kind of wideband reverse-time migration imaging method of marine adaptive compacting ghost reflection | |
| CN108363097B (en) | A kind of seismic data offset imaging method | |
| US8103453B2 (en) | Method of seismic data interpolation by projection on convex sets | |
| CN107678062B (en) | The integrated forecasting deconvolution of the domain hyperbolic Radon and feedback loop methodology multiple suppression model building method | |
| Kim et al. | Frequency-domain reverse-time migration with source estimation | |
| CN108108331B (en) | A kind of finite difference formulations method based on quasi- spatial domain equations for elastic waves | |
| CN108051855B (en) | A kind of finite difference formulations method based on quasi- spatial domain ACOUSTIC WAVE EQUATION | |
| EP3455650A1 (en) | Simultaneous source acquisition and separation on general related sampling grids | |
| CN110579795B (en) | Joint velocity inversion method based on passive source seismic waveform and reverse-time imaging thereof | |
| CN104459770B (en) | A kind of method for regularizing high-dimensional seismic data | |
| CN109541681A (en) | A kind of waveform inversion method of streamer seismic data and a small amount of OBS data aggregate | |
| CN110174702A (en) | A kind of method and system that marine seismic data low frequency weak signal is restored | |
| Yang et al. | An optimal implicit staggered-grid finite-difference scheme based on the modified Taylor-series expansion with minimax approximation method for elastic modeling | |
| CN105242313B (en) | A kind of bearing calibration of elastic wave reverse-time migration polarity inversion and system | |
| Xiang‐Bo et al. | Anisotropic Radon transform and its application to demultiple | |
| Gong et al. | Combined migration velocity model-building and its application in tunnel seismic prediction | |
| CN106842300A (en) | A kind of high efficiency multi-component seismic data true amplitude migration imaging method | |
| Qu et al. | 3-D Least-Squares Reverse Time Migration in Curvilinear-τ Domain | |
| WO2015155597A2 (en) | Attenuating pseudo s-waves in acoustic anisotropic wave propagation | |
| CN109557581A (en) | Reconstruction of seismic data method and system based on Fourier transformation | |
| Ying et al. | Denoise investigation on prestack reverse time migration based on GPU/CPU collaborative parallel accelerating computation | |
| CN108680957B (en) | Local cross-correlation time-frequency domain Phase-retrieval method based on weighting |
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 |