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
- longitudinal
- transverse
- transverse wave
- vsp
- 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 68
- 238000000926 separation method Methods 0.000 claims abstract description 23
- 238000005457 optimization Methods 0.000 claims abstract description 10
- 238000004088 simulation Methods 0.000 claims abstract description 9
- 238000001914 filtration Methods 0.000 claims abstract description 7
- 239000002245 particle Substances 0.000 claims description 12
- 230000005540 biological transmission Effects 0.000 claims description 9
- 238000012545 processing Methods 0.000 claims description 8
- 238000009499 grossing Methods 0.000 claims description 6
- 239000011159 matrix material Substances 0.000 claims description 6
- 230000008878 coupling Effects 0.000 claims description 4
- 238000010168 coupling process Methods 0.000 claims description 4
- 238000005859 coupling reaction Methods 0.000 claims description 4
- 230000010287 polarization Effects 0.000 claims description 4
- 238000004422 calculation algorithm Methods 0.000 claims description 3
- 238000001514 detection method Methods 0.000 claims description 3
- 230000002068 genetic effect Effects 0.000 claims description 3
- 238000006243 chemical reaction Methods 0.000 claims description 2
- 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 3
- 230000000694 effects Effects 0.000 description 3
- 238000003384 imaging method Methods 0.000 description 3
- 238000010276 construction Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 230000002159 abnormal effect Effects 0.000 description 1
- 238000010521 absorption reaction Methods 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000005422 blasting Methods 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000007598 dipping method Methods 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 230000003313 weakening effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/282—Application of seismic models, synthetic seismograms
-
- 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 the technical field of exploration geophysics, and relates to a novel full waveform inversion method which can be used for VSP seismic data inversion to obtain a high-precision converted wave velocity model.
Background
With the development of oil and gas exploration and the continuous improvement of seismic technology in wells, the research and application of VSP (vertical seismic profiling) exploration technology are rapidly developed. In areas of complex construction or complex media, VSP technology is effective in understanding the characteristics of geological structures, petrophysical properties, etc. near the wellbore. By combining VSP imaging with ground three-dimensional seismic data and lithological data, fine seismic attribute analysis can be performed, key positions can be further verified, contact relation between rock stratums can be cleared, an accurate geological model can be established, research on heterogeneity of the stratums can be performed, and technical support is provided for adjustment of development schemes. Compared with the ground earthquake, the VSP data contains more abundant wave field information, has high resolution and high signal-to-noise ratio, and the seismic record received in the well only passes through a low deceleration zone once, so that the attenuation of the seismic wave energy is smaller compared with the ground data, and the VSP data has the specific advantage of deep imaging. Secondly, due to observation limitation, the ground seismic exploration hardly observes the information of the target stratum salt hill body and the lateral wing and the salt. VSP technology can make up for these deficiencies in surface data. Therefore, VSP technology has its own advantages and effects, both in the exploration field and in the development field.
Velocity is a core parameter in seismic data processing flow. The high quality velocity modeling method plays a decisive role in subsequent offset imaging, lithology, lithofacies interpretation and reservoir prediction. So far, the data processing flow of VSP longitudinal wave data is relatively perfect, and the inversion accuracy of the longitudinal wave velocity is still good. However, due to the influence of many factors (such as extremely asymmetric transverse wave propagation, cross-coupled longitudinal and transverse waves, transverse wave anisotropy, etc.), the inversion accuracy of the converted wave velocity is still low. Elastic Full Waveform Inversion (EFWI) has the potential to obtain precise compressional, shear and density and is therefore becoming increasingly popular in exploration geophysics. However, due to multi-parameter coupling, the full waveform inversion of the elastic wave has serious nonlinearity and strong multi-solution property, and cannot be popularized in a large scale in the industry. In order to eliminate the influence of the longitudinal wave velocity on the inversion accuracy of the transverse wave velocity, it is necessary to further study the inversion method of the converted wave full waveform.
The invention provides a novel full waveform inversion method of VSP seismic data converted waves. Aiming at VSP observation modes (surface blasting and in-well receiving), PS wave (converted wave) records are separated through a global optimization parameter inversion method, and nonlinearity and multi-solution are relieved by adopting the separated converted waves for inversion; the method comprises the steps of obtaining longitudinal and transverse wave fields at different moments by solving an elastic wave equation and an adjoint equation of longitudinal and transverse wave separation, calculating the gradient of a target function relative to the transverse wave velocity by adopting a forward longitudinal wave field and an adjoint transverse wave field, and weakening the crosstalk effect among different wave modes to improve the inversion accuracy of the transverse wave velocity.
Disclosure of Invention
The invention aims to provide a full waveform inversion method for VSP seismic data converted waves.
The technical scheme adopted by the invention is carried out according to the following steps:
(1) separating PS waves (converted waves) of original VSP (vertical seismic profile) seismic data;
separating the uplink PS converted wave by a global optimization parameter inversion method;
(2) establishing a transverse wave velocity initial model;
inputting a longitudinal wave velocity model, obtaining a transverse wave velocity model according to an empirical formula, and smoothing the transverse wave velocity model to obtain an initial transverse wave velocity model;
(3) carrying out forward continuation on a seismic source wave field, and separating a longitudinal and transverse wave field;
and obtaining a positive transmission longitudinal and transverse wave field of each time step and a simulated VSP PS wave record by solving an elastic wave equation of longitudinal and transverse wave separation.
(4) Converting wave residual error reverse continuation, and separating longitudinal and transverse wave fields;
obtaining a reverse transmission longitudinal and transverse wave field of each time step by solving an adjoint elastic wave equation of longitudinal and transverse wave separation;
(5) and calculating the velocity gradient of the shear wave.
Calculating the shear wave velocity gradient by adopting the forward-transmitted longitudinal wave field and the backward-transmitted shear wave field obtained in the step (3) and the step (4);
(6) calculating an iteration step length;
solving an iteration step length by adopting a three-point parabola fitting method;
(7) updating the transverse wave velocity model;
when the convergence condition is not met, repeating the steps (3) to (6), and when the convergence condition is met, finishing the iteration;
(8) performing multi-scale inversion;
and performing frequency division processing on the VSP PS converted wave record by using wiener filtering, and gradually performing inversion from low frequency to high frequency.
Specifically, the original VSP seismic data PS converted wave separation method in the VSP seismic data converted wave full waveform inversion step (1) is as follows:
the invention discloses a parameter inversion method, which is a VSP data wave field separation method based on the apparent velocity and polarization characteristics of seismic waves.
Establishing a minimized target functional of observation data and simulation records for each depth point:
therein, ΨPdown,ΨPup,ΨSdownAnd ΨSupIncident angles v of down longitudinal wave, up longitudinal wave, down transverse wave and up transverse wave respectivelyPAnd vSIs the longitudinal and transverse wave velocity, u and uobsFor simulating and observing the wavefield, ω is the angular frequency, [ ω ]1,ω1]Is the frequency range.
And (3) obtaining the optimal propagation speed and propagation angle through the iteration of a global optimization genetic algorithm. And solving the equation (2) to realize the separation of the downlink longitudinal wave, the uplink longitudinal wave, the downlink transverse wave and the uplink transverse wave.
u(ω)＝[u(-M1,ω)…u(0,ω)…u(M2,ω)]T
＝aPdownuPdown(ω)+aPupuPup(ω)+aSdownuSdown(ω)+aSupuSup(ω), (2)
Wherein [ u ]Pdown,uPup,uSdown,uSup]Respectively, a down-going longitudinal wave, an up-going longitudinal wave, a down-going transverse wave and an up-going transverse wave, [ a ]Pdown,aPup,aSdown,aup]For the corresponding polarization vector, [ -M1,…,0,…,M2]The number of the detection point in the space window.
When the seismic source is a longitudinal wave source, the separated up-going transverse wave uSupThe inverse fourier transform (in the frequency domain) yields the observation PS converted wave record (decomposed into horizontal and vertical components according to the propagation angle).
The method for establishing the shear wave velocity initial model in the step (2) is as follows:
longitudinal and transverse wave velocity empirical formula:
and obtaining a transverse wave velocity model by adopting an empirical formula according to the input longitudinal wave velocity model, and smoothing the transverse wave velocity model through a Gaussian window to obtain an initial transverse wave velocity model.
In the step (3), the seismic source wave field is extended in the positive direction, and the method for separating the longitudinal wave field from the transverse wave field is as follows:
first order elastic wave velocity-stress equation:
wherein λ and μ are Lame constants, ρ is density, (v)x,vz) Is the particle vibration velocity vector (tau)xx,τzz,τxz) Is the stress vector.
The invention adopts the speed-stress elastic wave equation of longitudinal and transverse wave separation:
wherein (v)xp,vzp) Is the vibration velocity vector of longitudinal wave particle, (v)xs,vzs) Is the transverse wave particle vibration velocity vector (tau)xxp,τzzp) Is a longitudinal wave stress vector (tau)xxs,τzzs,τxzs) Is the shear wave stress vector.
And (3) solving the equation (6) by adopting a high-order staggered grid finite difference numerical value to obtain a positive transmission longitudinal and transverse wave field of each time step and simulating VSP PS wave record.
Step (4) reverse continuation of converted wave residual errors, and the method for separating longitudinal and transverse wave fields comprises the following steps:
in the adjoint state method, the adjoint operator satisfies:
<Lu,u*>＝<u,L*u*>, (7)
wherein: l is the positive operator, u is the particle velocity and stress vector, and "", is its accompanying state. Substituting equation (4) into equation (7) yields the conventional adjoint equation:
similarly, the adjoint equation of longitudinal and transverse wave separation can be obtained based on equation (6):
wherein ". sup." denotes the adjoint wavefield of the respective wavefield,for observation of the coupling of longitudinal and transverse waves and for simulation of VSP seismic recordings,converted-wave seismic records are isolated observations and simulated VSP.
And (3) taking the record residual errors of the simulated and observed VSP PS waves as a seismic source, and solving equations (10) and (11) by adopting a high-order staggered finite difference numerical value to obtain a reverse-propagation longitudinal and transverse wave field of each time step.
The method for calculating the shear wave velocity gradient in the step (5) comprises the following steps:
the shear wave velocity gradient in the conventional full waveform inversion is as follows:
wherein E is a target functional.
The invention adopts forward-transmitted longitudinal wave field and backward-transmitted transverse wave field to calculate the transverse wave velocity gradient:
wherein,
by adopting a new gradient formula (equation 13), the interference of the longitudinal wave velocity can be relieved, and the inversion accuracy of the transverse wave velocity is improved.
The iterative step length calculating method in the step (6) is as follows:
in the invention, a three-point parabolic interpolation method is adopted to calculate the step length. And selecting three heuristic points, and constructing a parabola by using the three points, wherein the extreme point of the parabola is the optimal step length.
Three heuristic step sizes a0,a1,a2Corresponding target functional is Es0,Es1,Es2Then the optimal step length is:
the method for updating the shear wave velocity model in the step (7) is as follows:
updating the transverse wave velocity by an L-BFGS method:
wherein m iskAnd mk+1Shear wave velocity models, H, for the current iteration and the next iteration, respectivelykAn approximation matrix that is the inverse of the hessian matrix.
The multi-scale inversion method in step 8 is as follows:
and performing frequency division processing on the VSP PS converted wave record by using wiener filtering, and gradually performing inversion from low frequency to high frequency. The wiener filter is:
wherein f iswienerIs a wiener filter, WoriginalIs a primary dominant frequency wavelet, WtargetThe target dominant frequency wavelet, ω is the angular frequency, and ε is the stationary coefficient, representing the complex conjugate.
And (4) taking the inversion result with large scale as an initial model of adjacent inversion with small scale, and repeating the steps (3) to (7) until a convergence condition is met.
The method has the beneficial effect that the full waveform inversion accuracy of the VSP seismic data converted wave is improved by a longitudinal and transverse wave separation strategy. The invention aims to make up for the defect of illumination limitation of a ground observation system, improve the modeling precision of the transverse wave velocity, and particularly improve the inversion effect of a salt dome flank and a steep dip angle target layer.
Drawings
FIG. 1 is a flow chart of a VSP seismic data converted wave full waveform inversion method;
FIG. 2 a complex construction model;
FIG. 3 shows a shear velocity gradient at a first iteration at a first scale;
FIG. 4 shows the inversion results of shear wave velocity;
FIG. 5 is a comparison graph of lateral longitudinal lateral transverse wave velocity curves of the well;
FIG. 6 transforms the shear wave target functional convergence curve.
Detailed Description
The present invention will be described in detail with reference to the following embodiments.
As shown in fig. 1, a flow chart for implementing full waveform inversion of VSP seismic data converted waves in the present invention specifically includes:
(1) separating PS converted waves of original VSP seismic data;
separating the uplink PS converted wave by a global optimization parameter inversion method;
(2) establishing a transverse wave velocity initial model;
inputting a longitudinal wave velocity model, obtaining a transverse wave velocity model according to an empirical formula, and smoothing the transverse wave velocity model to obtain an initial transverse wave velocity model;
(3) carrying out forward continuation on a seismic source wave field, and separating a longitudinal and transverse wave field;
and obtaining a positive transmission longitudinal and transverse wave field of each time step and a simulated VSP PS wave record by solving an elastic wave equation of longitudinal and transverse wave separation.
(4) Converting wave residual error reverse continuation, and separating longitudinal and transverse wave fields;
obtaining a reverse transmission longitudinal and transverse wave field of each time step by solving an adjoint elastic wave equation of longitudinal and transverse wave separation;
(5) and calculating the velocity gradient of the shear wave.
Calculating the shear wave velocity gradient by adopting the forward-transmitted longitudinal wave field and the backward-transmitted shear wave field obtained in the step (3) and the step (4);
(6) calculating an iteration step length;
solving an iteration step length by adopting a three-point parabola fitting method;
(7) updating the transverse wave velocity model;
when the convergence condition is not met, repeating the steps (3) to (6), and when the convergence condition is met, finishing the iteration;
(8) performing multi-scale inversion;
and performing frequency division processing on the VSP PS converted wave record by using wiener filtering, and gradually performing inversion from low frequency to high frequency.
The PS conversion wave separation method of the original VSP seismic data in the step (1) is as follows:
in the invention, a global optimization parameter inversion method is adopted to separate PS converted waves.
Establishing a minimized target functional of observation data and simulation records for each depth point:
therein, ΨPdown,ΨPup,ΨSdownAnd ΨSupIncident angles v of down longitudinal wave, up longitudinal wave, down transverse wave and up transverse wave respectivelyPAnd vSIs the longitudinal and transverse wave velocity, u and uobsFor simulating and observing the wavefield, ω is the angular frequency, [ ω ]1,ω1]Is the frequency range.
And (3) obtaining the optimal propagation speed and propagation angle through the iteration of a global optimization genetic algorithm. And solving the equation (2) to realize the separation of the downlink longitudinal wave, the uplink longitudinal wave, the downlink transverse wave and the uplink transverse wave.
u(ω)＝[u(-M1,ω)…u(0,ω)…u(M2,ω)]T
＝aPdownuPdown(ω)+aPupuPup(ω)+aSdownuSdown(ω)+aSupuSup(ω), (2)
Wherein [ u ]Pdown,uPup,uSdown,uSup]Respectively, a down-going longitudinal wave, an up-going longitudinal wave, a down-going transverse wave and an up-going transverse wave, [ a ]Pdown,aPup,aSdown,aup]For the corresponding polarization vector, [ -M1,…,0,…,M2]The number of the detection point in the space window.
When the seismic source is a longitudinal wave source, the separated up-going transverse wave uSupThe inverse fourier transform (in the frequency domain) yields the observation PS converted wave record (decomposed into horizontal and vertical components according to the propagation angle).
The method for establishing the shear wave velocity initial model in the step (2) is as follows:
longitudinal and transverse wave velocity empirical formula:
and obtaining a transverse wave velocity model by adopting an empirical formula according to the input longitudinal wave velocity model, and smoothing the transverse wave velocity model through a Gaussian window to obtain an initial transverse wave velocity model.
In the step (3), the seismic source wave field is extended in the positive direction, and the method for separating the longitudinal wave field from the transverse wave field is as follows:
first order elastic wave velocity-stress equation:
wherein λ and μ are Lame constants, ρ is density, (v)x,vz) Is the particle vibration velocity vector (tau)xx,τzz,τxz) Is the stress vector.
The invention adopts the speed-stress elastic wave equation of longitudinal and transverse wave separation:
wherein (v)xp,vzp) Is the vibration velocity vector of longitudinal wave particle, (v)xs,vzs) Is the transverse wave particle vibration velocity vector (tau)xxp,τzzp) Is a longitudinal wave stress vector (tau)xxs,τzzs,τxzs) Is the shear wave stress vector.
And (3) solving the equation (6) by adopting a high-order staggered grid finite difference numerical value to obtain a positive transmission longitudinal and transverse wave field of each time step and simulating VSP PS wave record.
Step (4) reverse continuation of converted wave residual errors, and the method for separating longitudinal and transverse wave fields comprises the following steps:
in the adjoint state method, the adjoint operator satisfies:
<Lu,u*>＝<u,L*u*>, (7)
wherein: l is the positive operator, u is the particle velocity and stress vector, and "", is its accompanying state. Substituting equation (4) into equation (7) yields the conventional adjoint equation:
similarly, the adjoint equation of longitudinal and transverse wave separation can be obtained based on equation (6):
wherein ". sup." denotes the adjoint wavefield of the respective wavefield,for observation of the coupling of longitudinal and transverse waves and for simulation of VSP seismic recordings,converted-wave seismic records are isolated observations and simulated VSP.
And (3) taking the record residual errors of the simulated and observed VSP PS waves as a seismic source, and solving equations (10) and (11) by adopting a high-order staggered finite difference numerical value to obtain a reverse-propagation longitudinal and transverse wave field of each time step.
The method for calculating the shear wave velocity gradient in the step (5) comprises the following steps:
the shear wave velocity gradient in the conventional full waveform inversion is as follows:
wherein E is a target functional.
The invention adopts forward-transmitted longitudinal wave field and backward-transmitted transverse wave field to calculate the transverse wave velocity gradient:
wherein,
by adopting a new gradient formula (equation 13), the interference of the longitudinal wave velocity can be relieved, and the inversion accuracy of the transverse wave velocity is improved.
The iterative step length calculating method in the step (6) is as follows:
in the invention, a three-point parabolic interpolation method is adopted to calculate the step length. And selecting three heuristic points, and constructing a parabola by using the three points, wherein the extreme point of the parabola is the optimal step length.
Three heuristic step sizes a0,a1,a2Corresponding target functional is Es0,Es1,Es2Then the optimal step length is:
the method for updating the shear wave velocity model in the step (7) is as follows:
updating the transverse wave velocity by an L-BFGS method:
wherein m iskAnd mk+1Shear wave velocity models, H, for the current iteration and the next iteration, respectivelykAn approximation matrix that is the inverse of the hessian matrix.
The multi-scale inversion method in the step (8) is as follows:
and performing frequency division processing on the VSP PS converted wave record by using wiener filtering, and gradually performing inversion from low frequency to high frequency. The wiener filter is:
wherein f iswienerIs a wiener filter, WoriginalIs a primary dominant frequency wavelet, WtargetThe target dominant frequency wavelet, ω is the angular frequency, and ε is the stationary coefficient, representing the complex conjugate.
And (4) taking the inversion result with large scale as an initial model of adjacent inversion with small scale, and repeating the steps (3) to (7) until a convergence condition is met.
Due to the adoption of the technical scheme, the invention has the following advantages: 1. the modeling precision of the shear wave velocity model can be improved. 2. The inversion accuracy of complex structures such as salt dome flanks and steeply dipping target layers can be improved.
The effectiveness of the inversion method of the full waveform of the converted wave of the VSP seismic data proposed in the present invention is analyzed by an example.
In fig. 2, (a) - (c) show the true compressional velocity, true shear velocity and initial shear velocity model of a complex structural model, respectively. The calculation area was 200 × 200 grid points, the spatial interval was 10m, the time step was 1ms, and the total recording time was 1.7 s. The source is a 15Hz Rake wavelet, which is placed under positive stress (longitudinal source). There is a high speed body at the lower right of the model, and the longitudinal wave velocity and the transverse wave velocity are different at the high speed body. The well is located at x 1980m and 180 receivers are placed in the well, one every 10 meters from z 100 m. The shot points are located on the earth surface, and 26 shots are taken every 70m from x-150 m. Numerical simulation is carried out by adopting a finite difference method of time second order and space tenth order precision, and boundary reflection is eliminated by a PML absorption boundary. And performing multi-scale inversion by adopting two scales, wherein the first scale is 15Hz, and the second scale is 30 Hz. FIG. 3 is a shear velocity gradient for a first iteration of a first scale. It can be seen that the shallow layer beside the well has clear structure, the deep layer has accurate depiction of the lateral wing of the salt dome. Fig. 4 (a) and (b) show inversion results of the first-scale and second-scale shear wave velocities, respectively. As can be seen from the figure, the inversion accuracy of the full waveform inversion method for VSP seismic data converted waves provided by the invention is higher, the structures beside a well and on the salt dome side are clear and visible, and high-speed abnormal bodies are better recovered. Fig. 5 is a graph comparing the true shear wave velocity, the initial velocity and the inversion result at a position where x is 1650 m. The inversion result is well matched with the real result. Fig. 6 shows a convergence curve of the objective function Es, which further shows the effectiveness of the VSP seismic data converted wave full waveform inversion method of the present invention.
The invention relates to a novel full waveform inversion method for VSP seismic data converted waves, which can obtain a high-precision transverse wave velocity model; the inversion accuracy of complex geological structures (such as salt dome flanks and steep dip interfaces) can be greatly improved.
The above description is only for the purpose of illustrating the preferred embodiments of the present invention and is not intended to limit the present invention in any way, and all simple modifications, equivalent variations and modifications made to the above embodiments according to the technical spirit of the present invention are within the scope of the present invention.
Claims (9)
1. A full waveform inversion method for VSP seismic data converted waves is characterized by comprising the following steps:
(1) separating PS converted waves of original VSP seismic data;
separating the uplink PS converted wave by a global optimization parameter inversion method;
(2) establishing a transverse wave velocity initial model;
inputting a longitudinal wave velocity model to obtain a transverse wave velocity model, and smoothing the transverse wave velocity model to obtain an initial transverse wave velocity model;
(3) carrying out forward continuation on a seismic source wave field, and separating a longitudinal and transverse wave field;
obtaining a positive transmission longitudinal and transverse wave field of each time step and a simulated VSPPS wave record by solving an elastic wave equation of longitudinal and transverse wave separation;
(4) converting wave residual error reverse continuation, and separating longitudinal and transverse wave fields;
obtaining a reverse transmission longitudinal and transverse wave field of each time step by solving an adjoint elastic wave equation of longitudinal and transverse wave separation;
(5) calculating the velocity gradient of the transverse wave;
calculating the shear wave velocity gradient by adopting the forward-transmission longitudinal wave field obtained in the step (3) and the backward-transmission shear wave field obtained in the step (4);
(6) calculating an iteration step length;
solving an iteration step length by adopting a three-point parabola fitting method;
(7) updating the transverse wave velocity model;
when the convergence condition is not met, repeating the steps (3) to (6), and when the convergence condition is met, finishing the iteration;
(8) performing multi-scale inversion;
and performing frequency division processing on the VSP PS converted wave record by using wiener filtering, and gradually performing inversion from low frequency to high frequency.
2. The VSP seismic data converted waveform full waveform inversion method of claim 1, wherein: the PS conversion wave separation method of the original VSP seismic data in the step (1) is as follows:
separating PS converted waves by adopting a global optimization parameter inversion method, and establishing a minimized target functional of observation data and simulation record for each depth point:
therein, ΨPdown,ΨPup,ΨSdownAnd ΨSupIncident angles v of down longitudinal wave, up longitudinal wave, down transverse wave and up transverse wave respectivelyPAnd vSIs the longitudinal and transverse wave velocity, u and uobsFor simulating and observing the wavefield, ω is the angular frequency, [ ω ]1,ω1]For the frequency range, the optimal propagation speed and propagation angle can be obtained through the iteration of a global optimization genetic algorithm; solving an equation (2) to realize separation of downlink longitudinal waves, uplink longitudinal waves, downlink transverse waves and uplink transverse waves;
u(ω)＝[u(-M1,ω)…u(0,ω)…u(M2,ω)]T
＝aPdownuPdown(ω)+aPupuPup(ω)+aSdownuSdown(ω)+aSupuSup(ω), (2)
wherein [ u ]Pdown,uPup,uSdown,uSup]Respectively, a down-going longitudinal wave, an up-going longitudinal wave, a down-going transverse wave and an up-going transverse wave, [ a ]Pdown,aPup,aSdown,aup]For the corresponding polarization vector, [ -M1,…,0,…,M2]The number of the detection point in the space window is the number, when the seismic source is a longitudinal wave source, the separated up-going transverse wave uSupAnd performing Fourier inversion to obtain an observation PS converted wave record.
3. The VSP seismic data converted waveform full waveform inversion method of claim 1, wherein: the method for establishing the shear wave velocity initial model in the step (2) is as follows:
longitudinal and transverse wave velocity empirical formula:
and obtaining a transverse wave velocity model by adopting an empirical formula according to the input longitudinal wave velocity model, and smoothing the transverse wave velocity model through a Gaussian window to obtain an initial transverse wave velocity model.
4. The VSP seismic data converted waveform full waveform inversion method of claim 1, wherein: the method for separating the longitudinal wave field from the transverse wave field in the step (3) comprises the following steps:
first order elastic wave velocity-stress equation:
wherein λ and μ are Lame constants, ρ is density, (v)x,vz) Is the particle vibration velocity vector (tau)xx,τzz,τxz) Is a stress vector;
adopting a speed-stress elastic wave equation of longitudinal and transverse wave separation:
wherein (v)xp,vzp) Is the vibration velocity vector of longitudinal wave particle, (v)xs,vzs) Is the transverse wave particle vibration velocity vector (tau)xxp,τzzp) Is a longitudinal wave stress vector (tau)xxs,τzzs,τxzs) Is the transverse wave stress vector;
and (3) solving the equation (6) by adopting a high-order staggered grid finite difference numerical value to obtain a positive transmission longitudinal and transverse wave field of each time step and simulating VSP PS wave record.
5. The VSP seismic data converted waveform full waveform inversion method of claim 1, wherein: the converted wave residual error reverse continuation in the step (4) and the method for separating the longitudinal wave field and the transverse wave field are as follows:
in the adjoint state method, the adjoint operator satisfies:
<Lu,u*>＝<u,L*u*>, (7)
wherein: l is the positive operator, u is the particle velocity and stress vector, "# is its companion state, substituting equation (4) into equation (7) yields the conventional companion equation:
obtaining a concomitant equation of longitudinal and transverse wave separation based on equation (6):
wherein, represents the adjoint wave field of the corresponding wave field,for observation of the coupling of longitudinal and transverse waves and for simulation of VSP seismic recordings,for separated observation and simulation VSP converted wave seismic records, simulating and observing VSP PS wave record residual errors are used as seismic sources, and equations (10) and (11) are solved by adopting high-order interleaved finite difference numerical values to obtain a reverse propagation longitudinal and transverse wave field of each time step.
6. The VSP seismic data converted waveform full waveform inversion method of claim 1, wherein: the method for calculating the shear wave velocity gradient in the step (5) comprises the following steps:
the shear wave velocity gradient in the conventional full waveform inversion is as follows:
wherein E is a target functional, and the shear wave velocity gradient is calculated by adopting a forward-transmission longitudinal wave field and a backward-transmission shear wave field:
wherein,
the new gradient formula is adopted to relieve the interference of longitudinal wave velocity, and further the inversion accuracy of transverse wave velocity is improved.
7. The VSP seismic data converted waveform full waveform inversion method of claim 1, wherein: the iterative step length calculating method in the step (6) is as follows:
selecting three trial points, and constructing a parabola by using the three points, wherein the extreme point of the parabola is the optimal step length, and the three trial step lengths a0,a1,a2Corresponding target functional is Es0,Es1,Es2Then the optimal step length is:
8. the VSP seismic data converted waveform full waveform inversion method of claim 1, wherein: the method for updating the shear wave velocity model in the step (7) is as follows:
updating the transverse wave velocity by an L-BFGS method:
wherein m iskAnd mk+1Shear wave velocity models, H, for the current iteration and the next iteration, respectivelykAn approximation matrix that is the inverse of the hessian matrix.
9. The VSP seismic data converted waveform full waveform inversion method of claim 1, wherein: the multi-scale inversion method in the step (8) is as follows:
performing frequency division processing on the VSP PS converted wave record by using wiener filtering, and gradually performing inversion from low frequency to high frequency, wherein the wiener filter is as follows:
wherein f iswienerIs a wiener filter, WoriginalIs a primary dominant frequency wavelet, WtargetAnd (3) repeating the steps (3) to (7) until a convergence condition is met, wherein omega is angular frequency, epsilon is a stable coefficient, represents complex conjugate, and the inversion result with a large scale is used as an initial model of adjacent small-scale inversion.
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 (11)
| 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 |
| CN109459789A (en) * | 2018-12-18 | 2019-03-12 | 吉林大学 | Time-domain full waveform inversion method based on amplitude decaying and linear interpolation |
| 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 |
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 |
|---|---|---|---|---|
| CN109459787B (en) * | 2018-10-09 | 2019-12-06 | 中国地质大学（武汉） | coal mine underground structure imaging method and system based on seismic channel wave full-waveform inversion |
| 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 |
| 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 | |
| CN108873066B (en) | Elastic medium wave equation reflected wave travel time inversion method | |
| CN105974470B (en) | A kind of multi-component seismic data least square reverse-time migration imaging method and system | |
| CN104597490B (en) | Multi-wave AVO reservoir elastic parameter inversion method based on accurate Zoeppritz equations | |
| CN103238158B (en) | Utilize the marine streamer data source inverting simultaneously that mutually related objects function is carried out | |
| CN102707317B (en) | Method of using absorption and attenuation characteristics of seismic wave for reservoir analysis | |
| CN101334483B (en) | Method for attenuating rayleigh wave scattered noise in earthquake data-handling | |
| CN109669212B (en) | Seismic data processing method, stratum quality factor estimation method and device | |
| CN102841376A (en) | Retrieval method for chromatography speed based on undulating surface | |
| CN111025387B (en) | Pre-stack earthquake multi-parameter inversion method for shale reservoir | |
| CN110187382B (en) | Traveling time inversion method for wave equation of reverse wave and reflected wave | |
| CN106842295A (en) | The waveform inversion method of logging information constrained | |
| CN110579795B (en) | Joint velocity inversion method based on passive source seismic waveform and reverse-time imaging thereof | |
| CN110261903B (en) | Underground seismic source passive positioning method based on reverse-time energy focusing | |
| CN107884829A (en) | A kind of method for combining compacting shallow sea OBC Multiple Attenuation in Seismic Data | |
| EP3167314B1 (en) | Method for obtaining estimates of a model parameter so as to characterise the evolution of a subsurface volume | |
| CN110737018B (en) | Method for modeling anisotropy of VSP seismic data | |
| CN104570095A (en) | Method for eliminating pseudo-reflection of inclined cable based on Radon conversion | |
| CN111983683A (en) | Prediction method and system for lake-facies limestone reservoir under low-well condition | |
| CN111025388A (en) | Multi-wave combined prestack waveform inversion method | |
| CN110873893A (en) | Deep water free surface multiple prediction and suppression method and system thereof | |
| CN113536638A (en) | High-precision seismic wave field simulation method based on discontinuous finite elements and staggered grids | |
| Ortiz-Aleman et al. | Imaging the Chicxulub central crater zone from large-scale seismic acoustic wave propagation and gravity modeling | |
| Yang et al. | Mitigating velocity errors in least-squares imaging using angle-dependent forward and adjoint Gaussian beam operators | |
| CN102466818A (en) | Method for imaging anisotropic medium by utilizing cross-well seismic data |
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 |