CN108845351A - A kind of VSP seismic data converted wave full waveform inversion method - Google Patents
A kind of VSP seismic data converted wave full waveform inversion method Download PDFInfo
- Publication number
- CN108845351A CN108845351A CN201810666468.2A CN201810666468A CN108845351A CN 108845351 A CN108845351 A CN 108845351A CN 201810666468 A CN201810666468 A CN 201810666468A CN 108845351 A CN108845351 A CN 108845351A
- Authority
- CN
- China
- Prior art keywords
- wave
- velocity
- shear wave
- longitudinal
- shear
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 70
- 238000000926 separation method Methods 0.000 claims abstract description 25
- 238000004088 simulation Methods 0.000 claims abstract description 13
- 238000005457 optimization Methods 0.000 claims abstract description 11
- 238000001914 filtration Methods 0.000 claims abstract description 7
- 239000002245 particle Substances 0.000 claims description 12
- 238000012545 processing Methods 0.000 claims description 7
- 239000011159 matrix material Substances 0.000 claims description 6
- 230000008878 coupling Effects 0.000 claims description 5
- 238000010168 coupling process Methods 0.000 claims description 5
- 238000005859 coupling reaction Methods 0.000 claims description 5
- 230000010287 polarization Effects 0.000 claims description 4
- 238000004422 calculation algorithm Methods 0.000 claims description 3
- 230000002068 genetic effect Effects 0.000 claims description 3
- 238000010304 firing Methods 0.000 claims 1
- 125000006850 spacer group Chemical group 0.000 claims 1
- 230000008569 process Effects 0.000 abstract description 3
- 230000009286 beneficial effect Effects 0.000 abstract description 2
- 150000003839 salts Chemical class 0.000 description 7
- 238000005516 engineering process Methods 0.000 description 5
- 238000011161 development Methods 0.000 description 4
- 230000008901 benefit Effects 0.000 description 3
- 238000006243 chemical reaction Methods 0.000 description 3
- 239000000463 material Substances 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 239000011435 rock Substances 0.000 description 2
- 208000035126 Facies Diseases 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 239000003989 dielectric material Substances 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 238000013316 zoning Methods 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
-
- 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/30—Analysis
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 invention discloses a kind of VSP seismic data converted wave full waveform inversion methods.Uplink PS converted wave is separated by global optimization parametric inversion method；Rule of thumb formula calculates shear wave velocity and is smoothed to obtain initial horizontal velocity model to it；Equations for elastic waves by solving longitudinal and shear wave separation obtains the main story of each time step wave field and simulation VSP PS wave record in length and breadth；Adjoint equations for elastic waves by solving longitudinal and shear wave separation obtains the anti-pass of each time step wave field in length and breadth；Using main story longitudinal wave field and anti-pass shear wave field computation shear wave velocity gradient；Iteration step length is sought using 3 curve-parabola-fitting methods；Shear wave velocity model is updated by L-BFGS method；Multi-scale inversion is carried out from low to high with Wiener filtering.The beneficial effects of the invention are as follows separating strategy by longitudinal and shear wave to improve the precision of VSP seismic data converted wave full waveform inversion, the efficiency of inverse process of complicated structure is improved.
Description
Technical field
The invention belongs to exploration geophysics technical fields, are related to a kind of full waveform inversion new method, with can be used for VSP
Data Inversion is shaken to obtain high-precision conversion velocity model.
Background technique
With the continuous improvement of seismic technology in the deep and well of oil-gas exploration and development, VSP (vertical seismic profiling (VSP)) exploration
The research and application of technology, have also obtained quick development.In complicated structure or complex dielectrics area, VSP technology can be effectively
Understand the features such as geologic structure, the physical properties of rock of near wellbore.Combined ground three dimensional seismic data and lithology money is imaged by VSP
Material, can be carried out fine seismic attributes analysis, can further verify key stratum position, get contact relation between rock stratum clear, builds
Accurate geological model is erected, the heterogeneity research on stratum is carried out, provides technical support for the adjustment of development plan.With ground
Seismic facies ratio, VSP data contain richer wave field information, have high resolution, signal-to-noise ratio is high, and the earthquake that downhole receiving arrives
Record is just through primary low deceleration strip, therefore the decaying for comparing surface data seismic wave energy is smaller, so having to imaging deep
Its peculiar advantage.Secondly surface seismic exploration is limited due to observation, be difficult to observe destination layer salt dome body and its flank and salt
Lower information.And VSP technology can make up these deficiencies of surface data.Therefore no matter in Exploration Domain, or in development field,
VSP technology has its distinctive advantage and effect.
Speed is the core parameter in seismic data process flow.The velocity modeling method of high quality is to subsequent inclined
Picture, lithology, petrographic interpretation and reservoir prediction is moved into play a decisive role.Up to the present, the data processing stream of VSP longitudinal wave data
Journey comparatively perfect, even can be with for the inversion accuracy of velocity of longitudinal wave.But (such as shear wave propagation pole is influenced due to receiving factors
Its asymmetry, longitudinal and shear wave coupling crosstalk, shear wave anisotropy etc.), converted wave velocity inversion precision is still lower.Elastic all-wave
Shape inverting (EFWI) has the potentiality for obtaining accurate velocity of longitudinal wave, shear wave velocity and density, therefore in exploration geophysics
Become to become more and more popular.But elastic wave full waveform inversion is due to Multi-parameter coupling, it is non-linear it is serious, multi-solution is strong, industrially
It could not also large-scale promotion.To eliminate influence of the velocity of longitudinal wave to S-wave velocity inversion precision, it is necessary to further research conversion
Wave full waveform inversion method.
The present invention provides a kind of VSP seismic data converted wave full waveform inversion new method.For VSP observed pattern (earth's surface
Blow out, downhole receiving), by global optimization parametric inversion method separation PS wave (converted wave) record, using the converted wave isolated
Inverting is carried out to alleviate non-linear and multi-solution；Difference is obtained by the equations for elastic waves and adjoint equation that solve longitudinal and shear wave separation
The wave field in length and breadth at moment, using positive longitudinal wave field and with ladder of the shear wave field computation objective function about shear wave velocity
Degree, weakens the crosstalk effect between different wave modes, to improve the inversion accuracy of shear wave velocity.
Summary of the invention
The purpose of the present invention is to provide a kind of VSP seismic data converted wave full waveform inversion methods.
The technical scheme adopted by the invention is that following the steps below：
(1) original VSP (vertical seismic profiling (VSP)) seismic data PS wave (converted wave) separation；
Uplink PS converted wave is separated by global optimization parametric inversion method；
(2) shear wave velocity initial model is established；
Velocity of longitudinal wave model is inputted, rule of thumb formula obtains shear wave velocity model, and is smoothed to obtain to it
Initial horizontal velocity model；
(3) source wavefield forward direction continuation, separating P wave from S wave wave field；
Equations for elastic waves by solving longitudinal and shear wave separation obtains main story wave field and the simulation in length and breadth of each time step
VSP PS wave record.
(4) converted wave residual error backward extension, separating P wave from S wave wave field；
Adjoint equations for elastic waves by solving longitudinal and shear wave separation obtains the anti-pass of each time step wave field in length and breadth；
(5) shear wave velocity gradient is calculated.
The main story longitudinal wave field and anti-pass shear wave field computation shear wave velocity gradient acquired using (3) step and (4) step；
(6) iteration step length is sought；
Iteration step length is sought using 3 curve-parabola-fitting methods；
(7) shear wave velocity model is updated；
When being unsatisfactory for the condition of convergence, repeats step (3) and arrive step (6), when meeting the condition of convergence, iteration terminates；
(8) Multi-scale inversion；
VSP PS converted wave is recorded with Wiener filtering and carries out scaling down processing, gradually carries out inverting from low to high.
Specifically, original VSP seismic data PS converts wavelength-division in VSP seismic data converted wave full waveform inversion step (1)
It is as follows from method：
Parametric inversion method is a kind of VSP data wave field separation method based on seismic wave apparent velocity and polarization characteristic, this hair
It is bright middle using global optimization parametric inversion method separation PS converted wave.
To each depth point, the minimization cost functional of observation data and analog record is established：
Wherein, ΨPdown,ΨPup,ΨSdownAnd ΨSupRespectively downlink longitudinal wave, uplink longitudinal wave, downlink shear wave and uplink are horizontal
The incidence angle of wave, vPAnd vSFor p-and s-wave velocity, u and uobsIt is angular frequency, [ω for simulation and observed wave field, ω1,ω1] it is frequency
Rate range.
Pass through global optimization genetic algorithm iteration available optimal spread speed and propagation angle.It solves equation (2)
Realize downlink longitudinal wave, uplink longitudinal wave, downlink shear wave and the separation of uplink shear wave.
U (ω)=[u (- M1,ω)…u(0,ω)…u(M2,ω)]T
=aPdownuPdown(ω)+aPupuPup(ω)+aSdownuSdown(ω)+aSupuSup(ω), (2)
Wherein, [uPdown,uPup,uSdown,uSup] it is respectively downlink longitudinal wave, uplink longitudinal wave, downlink shear wave and uplink shear wave
, [aPdown,aPup,aSdown,aup] it is corresponding polarization vector, [- M1,…,0,…,M2] be spatial window in geophone station number.
When focus is longitudinal wave source, the uplink shear wave u that isolatesSup(frequency domain) carries out Fourier inversion and just obtains
Observe PS converted wave record (horizontal component and vertical component are decomposed into according to propagation angle).
Shear wave velocity initial model method for building up is as follows in step (2)：
P-and s-wave velocity empirical equation：
According to the velocity of longitudinal wave model of input, shear wave velocity model is obtained using empirical equation, by Gaussian window to its into
Row smoothly obtains initial horizontal velocity model.
In step (3), source wavefield forward direction continuation, separating P wave from S wave wave field method is as follows：
Single order elastic wave velocity-stress equation：
Wherein, λ and μ is Lame constants, and ρ is density, (vx,vz) it is Particle Vibration Velocity vector, (τxx,τzz,τxz) it is to answer
Force vector.
Using speed-stress equations for elastic waves of longitudinal and shear wave separation in the present invention：
Wherein, (vxp,vzp) it is longitudinal-wave particle vibration velocity vector, (vxs,vzs) it is shear wave Particle Vibration Velocity vector,
(τxxp,τzzp) it is longitudinal wave stress vector, (τxxs,τzzs,τxzs) it is shear wave stress vector.
The main story of each time step wave field in length and breadth is obtained using high-order staggering mesh finite-difference numerical solution equation (6)
And simulation VSP PS wave record.
Step (4) converted wave residual error backward extension, separating P wave from S wave wave field method are as follows：
In Adjoint State Method, adjoint operator meets：
<Lu,u*>=<u,L*u*>, (7)
Wherein：L is forward operator, and u is particle velocity and stress vector, and " * " is it with state.Equation (4) is brought into
Equation (7) obtains conventional adjoint equation：
Similarly, the adjoint equation based on equation (6) available longitudinal and shear wave separation：
Wherein, " * " indicates the adjoint wave field of corresponding wave field,For longitudinal and shear wave coupling observation and
VSP earthquake record is simulated,For isolated observation and simulation VSP converted wave earthquake record.
To simulate and observe VSP PS wave record residual difference as focus, equation is solved using high-order limited difference numerical of interlocking
(10) and (11) obtain the anti-pass of each time step wave field in length and breadth.
It is as follows that step (5) calculates shear wave velocity gradient method：
Shear wave velocity gradient is in traditional full waveform inversion：
Wherein, E is cost functional.
The present invention uses main story longitudinal wave field and anti-pass shear wave field computation shear wave velocity gradient：
Wherein,
The interference of velocity of longitudinal wave can be alleviated by using new gradient formula (equation 13), and then improve shear wave velocity
Inversion accuracy.
It is as follows that step (6) seeks iteration step length method：
3 parabola interpolation method material calculations are used in the present invention.Three exploration points are chosen, are constructed with these three points
One parabola, this parabolical extreme point is both optimal step size.
Three exploration step-length a0,a1,a2Corresponding cost functional is Es0,Es1,Es2, then optimal step size be：
It is as follows that step (7) updates shear wave velocity model method：
Shear wave velocity is updated by L-BFGS method：
Wherein, mkAnd mk+1The respectively shear wave velocity model of current iteration and next iteration, HkIt is inverse for Hessian matrix
Approximate matrix.
Multi-scale inversion method is as follows in step 8：
VSP PS converted wave is recorded with Wiener filtering and carries out scaling down processing, gradually carries out inverting from low to high.
Wiener filter is：
Wherein, fwienerFor wiener filter, WoriginalIt is original dominant frequency wavelet, WtargetTarget dominant frequency wavelet, ω are
Angular frequency, ε are the coefficient of stability, and * indicates complex conjugate.
Initial model of the inversion result of large scale as adjacent small scale inversion repeats step (3) and arrives step (7), directly
To meeting the condition of convergence.
The beneficial effects of the invention are as follows separate strategy by longitudinal and shear wave to improve VSP seismic data converted wave full waveform inversion
Precision.The purpose of invention is the shortcomings that making up conventional surface observational system illumination limitation, to improve the modeling accuracy of shear wave velocity, especially
It is the efficiency of inverse process of salt dome flank and steep dip target zone.
Detailed description of the invention
Fig. 1 VSP seismic data converted wave full waveform inversion method flow chart；
Fig. 2 complicated structure model；
The shear wave velocity gradient of Fig. 3 the first scale first time iteration；
Fig. 4 S-wave velocity inversion result；
Longitudinal side line shear wave velocity curve comparison figure by Fig. 5 well；
Fig. 6 converted shear wave cost functional convergence curve.
Specific embodiment
The present invention is described in detail With reference to embodiment.
As shown in Figure 1, implementing the flow chart of VSP seismic data converted wave full waveform inversion for the present invention, specifically include：
(1) original VSP seismic data PS converted wave separation；
Uplink PS converted wave is separated by global optimization parametric inversion method；
(2) shear wave velocity initial model is established；
Velocity of longitudinal wave model is inputted, rule of thumb formula obtains shear wave velocity model, and is smoothed to obtain to it
Initial horizontal velocity model；
(3) source wavefield forward direction continuation, separating P wave from S wave wave field；
Equations for elastic waves by solving longitudinal and shear wave separation obtains main story wave field and the simulation in length and breadth of each time step
VSP PS wave record.
(4) converted wave residual error backward extension, separating P wave from S wave wave field；
Adjoint equations for elastic waves by solving longitudinal and shear wave separation obtains the anti-pass of each time step wave field in length and breadth；
(5) shear wave velocity gradient is calculated.
The main story longitudinal wave field and anti-pass shear wave field computation shear wave velocity gradient acquired using (3) step and (4) step；
(6) iteration step length is sought；
Iteration step length is sought using 3 curve-parabola-fitting methods；
(7) shear wave velocity model is updated；
When being unsatisfactory for the condition of convergence, repeats step (3) and arrive step (6), when meeting the condition of convergence, iteration terminates；
(8) Multi-scale inversion；
VSP PS converted wave is recorded with Wiener filtering and carries out scaling down processing, gradually carries out inverting from low to high.
Original VSP seismic data PS converted wave separation method is as follows in step (1)：
PS converted wave is separated using global optimization parametric inversion method in the present invention.
To each depth point, the minimization cost functional of observation data and analog record is established：
Wherein, ΨPdown,ΨPup,ΨSdownAnd ΨSupRespectively downlink longitudinal wave, uplink longitudinal wave, downlink shear wave and uplink are horizontal
The incidence angle of wave, vPAnd vSFor p-and s-wave velocity, u and uobsIt is angular frequency, [ω for simulation and observed wave field, ω1,ω1] it is frequency
Rate range.
Pass through global optimization genetic algorithm iteration available optimal spread speed and propagation angle.It solves equation (2)
Realize downlink longitudinal wave, uplink longitudinal wave, downlink shear wave and the separation of uplink shear wave.
U (ω)=[u (- M1,ω)…u(0,ω)…u(M2,ω)]T
=aPdownuPdown(ω)+aPupuPup(ω)+aSdownuSdown(ω)+aSupuSup(ω), (2)
Wherein, [uPdown,uPup,uSdown,uSup] it is respectively downlink longitudinal wave, uplink longitudinal wave, downlink shear wave and uplink shear wave
, [aPdown,aPup,aSdown,aup] it is corresponding polarization vector, [- M1,…,0,…,M2] be spatial window in geophone station number.
When focus is longitudinal wave source, the uplink shear wave u that isolatesSup(frequency domain) carries out Fourier inversion and just obtains
Observe PS converted wave record (horizontal component and vertical component are decomposed into according to propagation angle).
Shear wave velocity initial model method for building up is as follows in step (2)：
P-and s-wave velocity empirical equation：
According to the velocity of longitudinal wave model of input, shear wave velocity model is obtained using empirical equation, by Gaussian window to its into
Row smoothly obtains initial horizontal velocity model.
In step (3), source wavefield forward direction continuation, separating P wave from S wave wave field method is as follows：
Single order elastic wave velocity-stress equation：
Wherein, λ and μ is Lame constants, and ρ is density, (vx,vz) it is Particle Vibration Velocity vector, (τxx,τzz,τxz) it is to answer
Force vector.
Using speed-stress equations for elastic waves of longitudinal and shear wave separation in the present invention：
Wherein, (vxp,vzp) it is longitudinal-wave particle vibration velocity vector, (vxs,vzs) it is shear wave Particle Vibration Velocity vector,
(τxxp,τzzp) it is longitudinal wave stress vector, (τxxs,τzzs,τxzs) it is shear wave stress vector.
The main story of each time step wave field in length and breadth is obtained using high-order staggering mesh finite-difference numerical solution equation (6)
And simulation VSP PS wave record.
Step (4) converted wave residual error backward extension, separating P wave from S wave wave field method are as follows：
In Adjoint State Method, adjoint operator meets：
<Lu,u*>=<u,L*u*>, (7)
Wherein：L is forward operator, and u is particle velocity and stress vector, and " * " is it with state.Equation (4) is brought into
Equation (7) obtains conventional adjoint equation：
Similarly, the adjoint equation based on equation (6) available longitudinal and shear wave separation：
Wherein, " * " indicates the adjoint wave field of corresponding wave field,For longitudinal and shear wave coupling observation and
VSP earthquake record is simulated,For isolated observation and simulation VSP converted wave earthquake record.
To simulate and observe VSP PS wave record residual difference as focus, equation is solved using high-order limited difference numerical of interlocking
(10) and (11) obtain the anti-pass of each time step wave field in length and breadth.
It is as follows that step (5) calculates shear wave velocity gradient method：
Shear wave velocity gradient is in traditional full waveform inversion：
Wherein, E is cost functional.
The present invention uses main story longitudinal wave field and anti-pass shear wave field computation shear wave velocity gradient：
Wherein,
The interference of velocity of longitudinal wave can be alleviated by using new gradient formula (equation 13), and then improve shear wave velocity
Inversion accuracy.
It is as follows that step (6) seeks iteration step length method：
3 parabola interpolation method material calculations are used in the present invention.Three exploration points are chosen, are constructed with these three points
One parabola, this parabolical extreme point is both optimal step size.
Three exploration step-length a0,a1,a2Corresponding cost functional is Es0,Es1,Es2, then optimal step size be：
It is as follows that step (7) updates shear wave velocity model method：
Shear wave velocity is updated by L-BFGS method：
Wherein, mkAnd mk+1The respectively shear wave velocity model of current iteration and next iteration, HkIt is inverse for Hessian matrix
Approximate matrix.
Multi-scale inversion method is as follows in step (8)：
VSP PS converted wave is recorded with Wiener filtering and carries out scaling down processing, gradually carries out inverting from low to high.
Wiener filter is：
Wherein, fwienerFor wiener filter, WoriginalIt is original dominant frequency wavelet, WtargetTarget dominant frequency wavelet, ω are
Angular frequency, ε are the coefficient of stability, and * indicates complex conjugate.
Initial model of the inversion result of large scale as adjacent small scale inversion repeats step (3) and arrives step (7), directly
To meeting the condition of convergence.
The invention adopts the above technical scheme, which has the following advantages：1. building for shear wave velocity model can be improved
Mould precision.2. the inversion accuracy of complicated structure (such as salt dome flank and steep dip target zone) can be improved.
The VSP seismic data converted wave full waveform inversion method proposed in the present invention is analyzed below by an example
Validity.
(a)-(c) is respectively the true velocity of longitudinal wave of certain complicated structure model, true shear wave velocity and initial shear wave in Fig. 2
Rate pattern.Zoning is 200 × 200 mesh points, and space interval 10m, time step 1ms, the summary journal time is
1.7s.Focus is the Ricker wavelet of 15Hz, is added on direct stress (longitudinal wave source).There are a high speed body, and longitudinal wave in model lower right
Speed and the shear wave velocity form at the high speed body are different.Well location arranges 180 wave detectors, from z=at x=1980m in well
100m starts every 10 meters one.Shot point is located at earth's surface, every mono- big gun of 70m from x=150m, totally 26 big gun.Using time second order,
Ten rank precision finite difference calculus of space carries out numerical simulation, and PML absorbing boundary eliminates boundary reflection.It is carried out using two scales more
Scale inversion, the first scale are 15Hz, and the second scale is 30Hz.Fig. 3 is the shear wave velocity gradient of the first scale first time iteration.
It can be seen that the construction of shallow-layer is clear by well, infrastructure and salt dome flank are also accurately portrayed.(a), (b) in Fig. 4
The inversion result of respectively the first scale and the second scale shear wave velocity.As seen from the figure, VSP seismic data proposed by the present invention turns
The inversion accuracy for changing wave full waveform inversion method is higher, well take up salt dome flank construction it is high-visible, high-velocity anomaly body also obtains
Preferably restore.Fig. 5 is true shear wave velocity, initial velocity and inversion result contrast curve chart at x=1650m.Inversion result and
Legitimate reading preferably coincide.Fig. 6 provides the convergence curve of objective function Es, further demonstrates that VSP seismic data conversion of the present invention
The validity of wave full waveform inversion method.
The present invention is a kind of new VSP seismic data converted wave full waveform inversion method, can obtain high-precision shear wave
Rate pattern；The inversion accuracy of complex geological structure (such as salt dome flank and steep dip interface) can be greatly improved.
The above is only not to make limit in any form to the present invention to better embodiment of the invention
System, any simple modification that embodiment of above is made according to the technical essence of the invention, equivalent variations and modification,
Belong in the range of technical solution of the present invention.
Claims (9)
1. a kind of VSP seismic data converted wave full waveform inversion method, it is characterised in that follow the steps below：
(1) original VSP seismic data PS converted wave separation；
Uplink PS converted wave is separated by global optimization parametric inversion method；
(2) shear wave velocity initial model is established；
Velocity of longitudinal wave model is inputted, obtains shear wave velocity model, and be smoothed to obtain initial horizontal velocity model to it；
(3) source wavefield forward direction continuation, separating P wave from S wave wave field；
Equations for elastic waves by solving longitudinal and shear wave separation obtains the main story of each time step wave field and the VSP of simulation in length and breadth
PS wave record；
(4) converted wave residual error backward extension, separating P wave from S wave wave field；
Adjoint equations for elastic waves by solving longitudinal and shear wave separation obtains the anti-pass of each time step wave field in length and breadth；
(5) shear wave velocity gradient is calculated；
The anti-pass shear wave field computation shear wave velocity gradient that the main story longitudinal wave field and step (4) acquired using step (3) is acquired；
(6) iteration step length is sought；
Iteration step length is sought using 3 curve-parabola-fitting methods；
(7) shear wave velocity model is updated；
When being unsatisfactory for the condition of convergence, repeats step (3) and arrive step (6), when meeting the condition of convergence, iteration terminates；
(8) Multi-scale inversion；
VSP PS converted wave is recorded with Wiener filtering and carries out scaling down processing, gradually carries out inverting from low to high.
2. according to VSP seismic data converted wave full waveform inversion method described in claim 1, it is characterised in that：The step (1)
In original VSP seismic data PS converted wave separation method it is as follows：
PS converted wave is separated using global optimization parametric inversion method, to each depth point, establishes observation data and analog record
Minimization cost functional：
Wherein, ΨPdown,ΨPup,ΨSdownAnd ΨSupRespectively downlink longitudinal wave, uplink longitudinal wave, downlink shear wave and uplink shear wave enters
Firing angle, vPAnd vSFor p-and s-wave velocity, u and uobsIt is angular frequency, [ω for simulation and observed wave field, ω1,ω1] it is frequency range,
Pass through global optimization genetic algorithm iteration available optimal spread speed and propagation angle；It solves equation (2) and realizes downlink
Longitudinal wave, uplink longitudinal wave, downlink shear wave and the separation of uplink shear wave；
U (ω)=[u (- M1,ω)…u(0,ω)…u(M2,ω)]T
=aPdownuPdown(ω)+aPupuPup(ω)+aSdownuSdown(ω)+aSupuSup(ω), (2)
Wherein, [uPdown,uPup,uSdown,uSup] it is respectively downlink longitudinal wave, uplink longitudinal wave, downlink shear wave and uplink shear wave wave field,
[aPdown,aPup,aSdown,aup] it is corresponding polarization vector, [- M1,…,0,…,M2] be spatial window in geophone station number, when
When focus is longitudinal wave source, the uplink shear wave u that isolatesSupIt carries out Fourier inversion and has just obtained observation PS converted wave record.
3. according to VSP seismic data converted wave full waveform inversion method described in claim 1, it is characterised in that：The step (2)
Middle shear wave velocity initial model method for building up is as follows：
P-and s-wave velocity empirical equation：
According to the velocity of longitudinal wave model of input, shear wave velocity model is obtained using empirical equation, it is carried out by Gaussian window flat
Cunning obtains initial horizontal velocity model.
4. according to VSP seismic data converted wave full waveform inversion method described in claim 1, it is characterised in that：The step (3)
Middle source wavefield forward direction continuation, separating P wave from S wave wave field method are as follows：
Single order elastic wave velocity-stress equation：
Wherein, λ and μ is Lame constants, and ρ is density, (vx,vz) it is Particle Vibration Velocity vector, (τxx,τzz,τxz) it is to answer force vector
Amount；
Speed-stress the equations for elastic waves separated using longitudinal and shear wave：
Wherein, (vxp,vzp) it is longitudinal-wave particle vibration velocity vector, (vxs,vzs) it is shear wave Particle Vibration Velocity vector, (τxxp,
τzzp) it is longitudinal wave stress vector, (τxxs,τzzs,τxzs) it is shear wave stress vector；
Using high-order staggering mesh finite-difference numerical solution equation (6) obtain the main story of each time step in length and breadth wave field and
Simulate VSP PS wave record.
5. according to VSP seismic data converted wave full waveform inversion method described in claim 1, it is characterised in that：The step (4)
Middle converted wave residual error backward extension, separating P wave from S wave wave field method are as follows：
In Adjoint State Method, adjoint operator meets：
<Lu,u*>=<u,L*u*>, (7)
Wherein：L is forward operator, and u is particle velocity and stress vector, and " * " is it with state, brings equation (4) into equation
(7) conventional adjoint equation is obtained：
The adjoint equation of longitudinal and shear wave separation is obtained based on equation (6)：
Wherein, * indicates the adjoint wave field of corresponding wave field,For the observation and simulation of longitudinal and shear wave coupling
VSP earthquake record,For isolated observation and VSP converted wave earthquake record is simulated, to simulate and see
Survey VSP PS wave record residual difference is focus, is solved equation (10) and (11) using interlock limited difference numerical of high-order and is obtained per a period of time
The anti-pass of spacer step wave field in length and breadth.
6. according to VSP seismic data converted wave full waveform inversion method described in claim 1, it is characterised in that：The step (5)
Middle calculating shear wave velocity gradient method is as follows：
Shear wave velocity gradient is in traditional full waveform inversion：
Wherein, E is cost functional, using main story longitudinal wave field and anti-pass shear wave field computation shear wave velocity gradient：
Wherein,
Alleviate the interference of velocity of longitudinal wave by using new gradient formula, and then improves the inversion accuracy of shear wave velocity.
7. according to VSP seismic data converted wave full waveform inversion method described in claim 1, it is characterised in that：The step (6)
In to seek iteration step length method as follows：
Three exploration points to be chosen, put one parabola of building with these three, this parabolical extreme point had both been optimal step size, three
A exploration step-length a0,a1,a2Corresponding cost functional is Es0,Es1,Es2, then optimal step size be：
8. according to VSP seismic data converted wave full waveform inversion method described in claim 1, it is characterised in that：The step (7)
Middle update shear wave velocity model method is as follows：
Shear wave velocity is updated by L-BFGS method：
Wherein, mkAnd mk+1The respectively shear wave velocity model of current iteration and next iteration, HkFor the approximation that Hessian matrix is inverse
Matrix.
9. according to VSP seismic data converted wave full waveform inversion method described in claim 1, it is characterised in that：The step (8)
Middle Multi-scale inversion method is as follows：
VSP PS converted wave is recorded with Wiener filtering and carries out scaling down processing, gradually carries out inverting, wiener from low to high
Filter is：
Wherein, fwienerFor wiener filter, WoriginalIt is original dominant frequency wavelet, WtargetTarget dominant frequency wavelet, ω are angular frequency
Rate, ε are the coefficient of stability, and * indicates complex conjugate, and initial model of the inversion result of large scale as adjacent small scale inversion repeats
Step (3) arrives step (7), until meeting the condition of convergence.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810666468.2A CN108845351A (en) | 2018-06-26 | 2018-06-26 | A kind of VSP seismic data converted wave full waveform inversion method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810666468.2A CN108845351A (en) | 2018-06-26 | 2018-06-26 | A kind of VSP seismic data converted wave full waveform inversion method |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN108845351A true CN108845351A (en) | 2018-11-20 |
Family
ID=64201961
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810666468.2A Pending CN108845351A (en) | 2018-06-26 | 2018-06-26 | A kind of VSP seismic data converted wave full waveform inversion method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108845351A (en) |
Cited By (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109459789A (en) * | 2018-12-18 | 2019-03-12 | 吉林大学 | Time-domain full waveform inversion method based on amplitude decaying and linear interpolation |
| CN109459787A (en) * | 2018-10-09 | 2019-03-12 | 中国地质大学（武汉） | Underground coal mine structure imaging method and system based on SEISMIC CHANNEL WAVE METHOD full waveform inversion |
| CN111751876A (en) * | 2020-07-07 | 2020-10-09 | 中油奥博(成都)科技有限公司 | Method and device for changing offset VSP (vertical seismic profiling) converted shear wave single-pass wave prestack depth offset |
| EP3761067A1 (en) * | 2019-07-01 | 2021-01-06 | CGG Services SAS | Full waveform inversion approach to building an s-wave velocity model using ps data |
| CN112773396A (en) * | 2021-01-13 | 2021-05-11 | 佟小龙 | Medical imaging method based on full waveform inversion, computer equipment and storage medium |
| CN113156498A (en) * | 2021-02-26 | 2021-07-23 | 中海石油(中国)有限公司 | Pre-stack AVO three-parameter inversion method and system based on homotopy continuation |
| CN113376695A (en) * | 2021-06-11 | 2021-09-10 | 中国矿业大学 | Full waveform inversion method suitable for complex collapse column of coal seam floor |
| CN113589366A (en) * | 2020-04-30 | 2021-11-02 | 中国石油化工股份有限公司 | Broadband fusion modeling method based on full waveform inversion |
| CN113987888A (en) * | 2021-11-16 | 2022-01-28 | 华中科技大学 | Bedrock incident-downlink wave number value separation method based on array observation |
| CN114660659A (en) * | 2022-03-30 | 2022-06-24 | 中国矿业大学 | Longitudinal and transverse wave decoupling high-precision microseism positioning method |
| CN116540297A (en) * | 2023-05-06 | 2023-08-04 | 中国科学院地质与地球物理研究所 | Full waveform inversion method, system and equipment for elastic wave seismic data |
| CN113987888B (en) * | 2021-11-16 | 2024-07-09 | 华中科技大学 | Bedrock incidence-downlink wave numerical value separation method based on array observation |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103499835A (en) * | 2013-10-13 | 2014-01-08 | 中国石油集团西北地质研究所 | Method for inverting near-surface velocity model by utilizing preliminary waveforms |
| CN106842295A (en) * | 2015-12-04 | 2017-06-13 | 中国石油化工股份有限公司 | The waveform inversion method of logging information constrained |
| WO2018035370A1 (en) * | 2016-08-19 | 2018-02-22 | Halliburton Energy Services, Inc. | Full waveform inversion of vertical seismic profile data for anisotropic velocities using pseudo-acoustic wave equations |
-
2018
- 2018-06-26 CN CN201810666468.2A patent/CN108845351A/en active Pending
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103499835A (en) * | 2013-10-13 | 2014-01-08 | 中国石油集团西北地质研究所 | Method for inverting near-surface velocity model by utilizing preliminary waveforms |
| CN106842295A (en) * | 2015-12-04 | 2017-06-13 | 中国石油化工股份有限公司 | The waveform inversion method of logging information constrained |
| WO2018035370A1 (en) * | 2016-08-19 | 2018-02-22 | Halliburton Energy Services, Inc. | Full waveform inversion of vertical seismic profile data for anisotropic velocities using pseudo-acoustic wave equations |
Non-Patent Citations (6)
| Title |
|---|
| LEONARDO GO ´MEZ ET AL.: "Wavefield separation methods for gradient filtering in time-domain full-waveform inversion", 《FIFTEENTH INTERNATIONAL CONGRESS OF THE BRAZILIAN GEOPHYSICAL SOCIETY》 * |
| MIN LOU ET AL.: "An improved parametric inversion methodology to separate P and Sv wavefields from VSP data", 《SEG HOUSTON 2013 ANNUAL MEETING》 * |
| ZHIMING REN ET AL.: "A hierarchical elastic full-waveform inversion scheme based on wavefield separation and the multistep-length approach", 《GEOPHYSICS》 * |
| 任志明: "声波和弹性波波动方程有限差分正反演方法研究", 《中国博士学位论文全文数据库 基础科学辑》 * |
| 孙思宇: "时间域弹性波全波形反演", 《中国优秀硕士学位论文全文数据库 基础科学辑》 * |
| 郭振波: "弹性介质波形反演方法研究", 《中国博士学位论文全文数据库 基础科学辑》 * |
Cited By (20)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109459787A (en) * | 2018-10-09 | 2019-03-12 | 中国地质大学（武汉） | Underground coal mine structure imaging method and system based on SEISMIC CHANNEL WAVE METHOD full waveform inversion |
| CN109459787B (en) * | 2018-10-09 | 2019-12-06 | 中国地质大学（武汉） | coal mine underground structure imaging method and system based on seismic channel wave full-waveform inversion |
| CN109459789A (en) * | 2018-12-18 | 2019-03-12 | 吉林大学 | Time-domain full waveform inversion method based on amplitude decaying and linear interpolation |
| EP3761067A1 (en) * | 2019-07-01 | 2021-01-06 | CGG Services SAS | Full waveform inversion approach to building an s-wave velocity model using ps data |
| US11215720B2 (en) | 2019-07-01 | 2022-01-04 | Cgg Services Sas | Full waveform inversion approach to building an S-wave velocity model using PS data |
| CN113589366A (en) * | 2020-04-30 | 2021-11-02 | 中国石油化工股份有限公司 | Broadband fusion modeling method based on full waveform inversion |
| CN113589366B (en) * | 2020-04-30 | 2023-10-20 | 中国石油化工股份有限公司 | Broadband fusion modeling method based on full waveform inversion |
| CN111751876A (en) * | 2020-07-07 | 2020-10-09 | 中油奥博(成都)科技有限公司 | Method and device for changing offset VSP (vertical seismic profiling) converted shear wave single-pass wave prestack depth offset |
| CN111751876B (en) * | 2020-07-07 | 2022-05-20 | 中油奥博(成都)科技有限公司 | Method and device for shifting prestack depth of converted shear wave single-pass wave of variable offset VSP (vertical seismic profiling) |
| CN112773396A (en) * | 2021-01-13 | 2021-05-11 | 佟小龙 | Medical imaging method based on full waveform inversion, computer equipment and storage medium |
| CN112773396B (en) * | 2021-01-13 | 2023-06-16 | 佟小龙 | Medical imaging method based on full waveform inversion, computer equipment and storage medium |
| CN113156498A (en) * | 2021-02-26 | 2021-07-23 | 中海石油(中国)有限公司 | Pre-stack AVO three-parameter inversion method and system based on homotopy continuation |
| CN113156498B (en) * | 2021-02-26 | 2024-01-26 | 中海石油(中国)有限公司 | Pre-stack AVO three-parameter inversion method and system based on homotopy continuation |
| CN113376695A (en) * | 2021-06-11 | 2021-09-10 | 中国矿业大学 | Full waveform inversion method suitable for complex collapse column of coal seam floor |
| CN113376695B (en) * | 2021-06-11 | 2022-07-05 | 中国矿业大学 | Full waveform inversion method suitable for complex collapse column of coal seam floor |
| CN113987888A (en) * | 2021-11-16 | 2022-01-28 | 华中科技大学 | Bedrock incident-downlink wave number value separation method based on array observation |
| CN113987888B (en) * | 2021-11-16 | 2024-07-09 | 华中科技大学 | Bedrock incidence-downlink wave numerical value separation method based on array observation |
| CN114660659A (en) * | 2022-03-30 | 2022-06-24 | 中国矿业大学 | Longitudinal and transverse wave decoupling high-precision microseism positioning method |
| CN116540297A (en) * | 2023-05-06 | 2023-08-04 | 中国科学院地质与地球物理研究所 | Full waveform inversion method, system and equipment for elastic wave seismic data |
| CN116540297B (en) * | 2023-05-06 | 2024-03-08 | 中国科学院地质与地球物理研究所 | Full waveform inversion method, system and equipment for elastic wave seismic data |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108845351A (en) | A kind of VSP seismic data converted wave full waveform inversion method | |
| CN107526101B (en) | A kind of acquisition and processing method obtaining earthquake reflected wave | |
| CN101334483B (en) | Method for attenuating rayleigh wave scattered noise in earthquake data-handling | |
| CN108802812A (en) | A kind of formation lithology inversion method of well shake fusion | |
| CN102645670B (en) | Observation system optimization design method based on stack response analysis | |
| CN101598811B (en) | Method for calculating shot-static correction for two-dimensional vertical seismic section data | |
| CN105334535A (en) | Method for identifying thin-reservoir subtle lithologic oil and gas reservoir | |
| CN102053270A (en) | Sedimentary formation unit-based seismic facies analysis method | |
| CN108873066A (en) | Elastic fluid fluctuates equation back wave Travel Time Inversion method | |
| CN106405651A (en) | Logging-matching-based method for constructing full-waveform inversion initial model | |
| CN106842295A (en) | The waveform inversion method of logging information constrained | |
| CN105277978A (en) | Method and device for determining near-ground-surface speed model | |
| CN101201409B (en) | Method for revising earthquake data phase | |
| CN107783187A (en) | It is a kind of that logging speed and seismic velocity are combined to the method for establishing 3D velocity field | |
| CN110187382B (en) | Traveling time inversion method for wave equation of reverse wave and reflected wave | |
| CN109669212A (en) | Seismic data processing technique, interval quality factors evaluation method and device | |
| CN107462924A (en) | A kind of absolute wave impedance inversion method independent of well-log information | |
| CN110007340A (en) | Salt dome speed density estimation method based on the direct envelope inverting of angle domain | |
| CN107728205B (en) | A kind of Formation pressure prediction method | |
| CN105607119B (en) | Near-surface model construction method and static correction value acquiring method | |
| CN107817516A (en) | Near surface modeling method and system based on preliminary wave information | |
| CN111983683A (en) | Prediction method and system for lake-facies limestone reservoir under low-well condition | |
| CN102053269A (en) | Analysis method of speed in seismic data | |
| CN112230274B (en) | While-drilling-oriented acoustic wave equation frequency domain reverse-time migration rapid imaging method | |
| Deng et al. | Hydrocarbon accumulation conditions and key exploration and development technologies for PL 19–3 oilfield |
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 | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20181120 |
|
| RJ01 | Rejection of invention patent application after publication |