CN107092029A - A kind of seismic inversion method and device - Google Patents
A kind of seismic inversion method and device Download PDFInfo
- Publication number
- CN107092029A CN107092029A CN201710280917.5A CN201710280917A CN107092029A CN 107092029 A CN107092029 A CN 107092029A CN 201710280917 A CN201710280917 A CN 201710280917A CN 107092029 A CN107092029 A CN 107092029A
- Authority
- CN
- China
- Prior art keywords
- reservoir
- vertical
- velocity
- wave
- mrow
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 59
- 230000008859 change Effects 0.000 claims abstract description 54
- 238000004364 calculation method Methods 0.000 claims description 8
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 2
- 239000000047 product Substances 0.000 description 20
- 238000001615 p wave Methods 0.000 description 7
- 230000008569 process Effects 0.000 description 6
- 239000002245 particle Substances 0.000 description 5
- 238000004088 simulation Methods 0.000 description 5
- 230000006854 communication Effects 0.000 description 4
- 239000011159 matrix material Substances 0.000 description 4
- 239000013598 vector Substances 0.000 description 4
- 230000005540 biological transmission Effects 0.000 description 3
- 238000004891 communication Methods 0.000 description 3
- 238000010168 coupling process Methods 0.000 description 3
- 238000005859 coupling reaction Methods 0.000 description 3
- 238000009795 derivation Methods 0.000 description 3
- 241000208340 Araliaceae Species 0.000 description 2
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 2
- 235000003140 Panax quinquefolius Nutrition 0.000 description 2
- 230000015572 biosynthetic process Effects 0.000 description 2
- 238000004422 calculation algorithm Methods 0.000 description 2
- 230000008878 coupling Effects 0.000 description 2
- 235000008434 ginseng Nutrition 0.000 description 2
- 238000011017 operating method Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- BVKZGUZCCUSVTD-UHFFFAOYSA-L Carbonate Chemical compound [O-]C([O-])=O BVKZGUZCCUSVTD-UHFFFAOYSA-L 0.000 description 1
- 230000003466 anti-cipated effect Effects 0.000 description 1
- 239000007795 chemical reaction product Substances 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000007405 data analysis Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000005315 distribution function Methods 0.000 description 1
- 238000005553 drilling Methods 0.000 description 1
- 235000013399 edible fruits Nutrition 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 238000007619 statistical method 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
-
- 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
- G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
- G01V2210/624—Reservoir parameters
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 provides a kind of seismic inversion method and device, wherein, this method includes：Based on velocity of longitudinal wave in reservoir to be measured when the value of vertical direction, velocity of longitudinal wave are in nearly vertical incidence with the situation of change at phase angle, horizontal slowness, velocity of longitudinal wave in horizontal and vertical difference value, determine the vertical slowness relational expression of reservoir plane compressional wave to be measured；Further according to shear wave velocity in reservoir to be measured in the value of vertical direction, the vertical slowness relational expression of reservoir transverse wave to be measured is determined；And the reflectance factor calculating mode of the new Method in Transverse Isotropic Medium with vertical symmetry axis is built, inverting is carried out to reservoir to be measured.In embodiments of the present invention, above-mentioned two formula not only has higher computational accuracy, also there is relatively simple form and relatively weak non-linear, reflectance factor after simplification calculates mode and higher computational accuracy is respectively provided with when strong anisotropy, strong impedance contrast and big incidence angle, and the surveying accuracy for carrying out oil field prospecting has obtained effective raising.
Description
Technical field
The present invention relates to technical field of geological exploration, more particularly to a kind of seismic inversion method and device.
Background technology
Seismic inversion also refers to the seismic data observed using earth's surface, with known geologic rule and drilling well, well logging
Data is constraint, thus it is speculated that go out the process of formation lithology construction.Divide from the seismic data used in seismic inversion, seismic inversion can divide
For：Prestack inversion and post-stack inversion；From inverting using the information of earthquake come point：Seismic inversion can be divided into：When seismic wave is travelled
Inverting and seismic amplitude inverting；Divide from the geology result of inverting, seismic inversion can be divided into：Construct inverting, wave impedance inversion,
Reservoir parameter inversion, Geostatistical Inversion etc..
At present, can typically be realized using following methods has the transverse isotropy of vertical symmetry axis in Seismic Reservoir
The inverting of (Vertical Transverse Isotropy, referred to as VTI) medium parameter：
1) inversion method based on the assumption of isotropy：
However, this method is unable to the VTI parameters that inverting obtains the shale reservoir with anisotropic character, and due to
Under VTI Anisotropic Conditions, isotropic reflectance factor calculating formula calculates obtained reflectance factor and there is larger error, makes
The common elastic parameters precision that must be obtained based on the inversion method inverting of the assumption of isotropy is not high.
2) full waveform inversion method：
Although this method can be using all-wave information come the VTI parameters of predicting reservoir.However, this method amount of calculation is huge
Greatly, the requirement of actual oil reservoir can not be met on inverting yardstick and computational efficiency.
3) amplitude based on VTI media with angle change (Amplitude Versus Angle, referred to as AVA) inverting
Method：
This method is compared with inversion method and full waveform inversion method based on the assumption of isotropy, and geological data is explained
Confidence level be improved, mainly can by combine VTI media accurate reflection coefficient formula and prestack trace gather on shake
Width the information such as changes with offset distance, to determine the VTI parameters in Seismic Reservoir (elastic parameter and anisotropic parameters).However,
During carrying out inverting, because accurate reflection coefficient approximate formula is extremely complex, it is difficult to be directly used in inverting.Meanwhile, by
It is relatively low in the accurate reflection coefficient formula computational accuracy of VTI media, in feelings such as strong anisotropy, strong impedance contrast and big incidence angles
Under condition, there is larger error in the VTI parameters that inverting is obtained, it is impossible to meet the fine sign requirement of actual oil reservoir.
The content of the invention
The invention provides a kind of seismic inversion method and device, to reach the inversion accuracy for improving VTI parameters, oil is instructed
The purpose of field Efficient Development.
The embodiments of the invention provide a kind of seismic inversion method, it can include：Existed based on velocity of longitudinal wave in reservoir to be measured
The value of vertical direction, the velocity of longitudinal wave are in nearly vertical incidence with the situation of change at phase angle, horizontal slowness, the compressional wave
Speed determines the vertical slowness relational expression of the reservoir plane compressional wave to be measured in horizontal and vertical difference value；According to storage to be measured
Shear wave velocity is closely hanging down in the value of vertical direction, the velocity of longitudinal wave in the value of vertical direction, the velocity of longitudinal wave in layer
With the situation of change at phase angle, horizontal slowness, the velocity of longitudinal wave in horizontal and vertical difference value when straight incident, it is determined that described
The vertical slowness relational expression of reservoir transverse wave to be measured；According to the vertical slowness relational expression of the reservoir plane compressional wave to be measured, institute
The vertical slowness relational expression of reservoir transverse wave to be measured is stated, the new Method in Transverse Isotropic Medium with vertical symmetry axis is built
Reflectance factor calculates mode；Mode is calculated using the reflectance factor, row amplitude is entered with angle change to the reservoir to be measured
Inverting.
In one embodiment, the vertical slowness of the reservoir plane compressional wave to be measured is determined, can be included：Treated described in calculating
Value of the velocity of longitudinal wave described in reservoir in vertical direction is surveyed, and calculates square reciprocal of the value, the first reservoir is used as
Parameter；Calculate square of the horizontal slowness, and by the horizontal slowness square with the velocity of longitudinal wave in nearly vertical incidence
When with phase angle situation of change be multiplied, be used as the second reservoir parameter；The velocity of longitudinal wave is calculated in horizontal and vertical difference
Value, the velocity of longitudinal wave in nearly vertical incidence with phase angle situation of change difference, as speed difference value, calculate described vertical
Wave velocity calculates the biquadratic of the horizontal slowness in square of the value of vertical direction, calculates the speed difference value, described
Value square, the difference between the biquadratic of the horizontal slowness, and regard resulting difference as the 3rd reservoir parameter；Meter
Calculate the root of the 3rd reservoir parameter, and using the root as the reservoir plane compressional wave to be measured vertical slowness.
In one embodiment, the vertical slowness of the reservoir plane compressional wave to be measured is determined according to below equation：
In above formula, qαRepresent the vertical slowness of the reservoir plane compressional wave to be measured, vP0Represent the velocity of longitudinal wave vertical
The value in direction, δ represents velocity of longitudinal wave situation of change with phase angle in nearly vertical incidence, and p represents that the level is slow
Degree, ε represents the velocity of longitudinal wave in horizontal and vertical difference value.
In one embodiment, the vertical slowness of the reservoir transverse wave to be measured is determined, can be included：Treated described in calculating
Survey shear wave velocity in reservoir and, in the inverse of the value of vertical direction, and calculate square reciprocal, be used as the 4th reservoir parameter；
Product square the biquadratic of the horizontal slowness between of the velocity of longitudinal wave in the value of vertical direction is calculated, is obtained
Five reservoir parameters；Calculate the velocity of longitudinal wave vertical direction value and the shear wave velocity the value of vertical direction ratio
Value square, and calculate the ratio square with the horizontal slowness square product as level value, by the ratio
Square be multiplied with the level value, obtain the 6th reservoir parameter；Calculate the 5th reservoir parameter and the 6th reservoir ginseng
Several differences, is used as the 7th reservoir parameter；Calculate the product of the speed difference value and the 7th reservoir parameter, and by gained
The result of product arrived is added with the 4th reservoir parameter, and the result after will add up is used as the 8th reservoir parameter；Calculate described
Eight reservoir parameters, the horizontal slowness square between difference, be used as the reservoir plane to be measured horizontal resulting difference
The vertical slowness of ripple.
In one embodiment, the vertical slowness of the reservoir transverse wave to be measured is determined according to below equation：
In above formula, qβRepresent the vertical slowness of the reservoir transverse wave to be measured, vS0Represent the shear wave velocity vertical
The value in direction, vP0The value of velocity of longitudinal wave in the reservoir to be measured in vertical direction is represented, p represents the horizontal slowness, ε tables
Show the velocity of longitudinal wave in horizontal and vertical difference value, δ represents the longitudinal wave velocity degree in nearly vertical incidence with phase angle
Situation of change.
The embodiment of the present invention additionally provides a kind of seismic inversion device, can include：Compressional wave determining module, for based on treating
Survey the change feelings with phase angle when the value of vertical direction, the velocity of longitudinal wave are in nearly vertical incidence of velocity of longitudinal wave in reservoir
Condition, horizontal slowness, the velocity of longitudinal wave determine the vertical slow of the reservoir plane compressional wave to be measured in horizontal and vertical difference value
Spend relational expression；Shear wave determining module, for the value according to shear wave velocity in reservoir to be measured in vertical direction, the velocity of longitudinal wave
Value, the velocity of longitudinal wave in vertical direction are in nearly vertical incidence with the situation of change at phase angle, horizontal slowness, described vertical
Wave velocity determines the vertical slowness relational expression of the reservoir transverse wave to be measured in horizontal and vertical difference value；Reflectance factor
Computing module, for shear wave in the vertical slowness relational expression according to the reservoir plane compressional wave to be measured, the reservoir plane to be measured
Vertical slowness relational expression, build the reflectance factor calculating mode of the new Method in Transverse Isotropic Medium with vertical symmetry axis；
Inverting module, for calculating mode using the reflectance factor, inverting of the row amplitude with angle change is entered to the reservoir to be measured.
In one embodiment, the compressional wave determining module can include：First reservoir parameter computing unit, for calculating
Velocity of longitudinal wave described in the reservoir to be measured and calculates square reciprocal of the value in the value of vertical direction, is used as
One reservoir parameter；Second reservoir parameter computing unit, square for calculating the horizontal slowness, and by the horizontal slowness
Square it is multiplied with the velocity of longitudinal wave in nearly vertical incidence with the situation of change at phase angle, is used as the second reservoir parameter；3rd
Reservoir parameter computing unit, is closely hanging down for calculating the velocity of longitudinal wave in horizontal and vertical difference value, the velocity of longitudinal wave
When straight incident with phase angle situation of change difference, be used as speed difference value, calculate velocity of longitudinal wave the taking in vertical direction
Square of value, calculates the biquadratic of the horizontal slowness, calculate the speed difference value, the value square, the level
Difference between the biquadratic of slowness, and it regard resulting difference as the 3rd reservoir parameter；The vertical slowness of compressional wave calculates single
Member, the root for calculating the 3rd reservoir parameter, and it regard the root as the vertical of the reservoir plane compressional wave to be measured
Slowness.
In one embodiment, the compressional wave determining module specifically can be used for determining the storage to be measured according to below equation
The vertical slowness of layer plane compressional wave：
In above formula, qαRepresent the vertical slowness of the reservoir plane compressional wave to be measured, vP0Represent the velocity of longitudinal wave vertical
The value in direction, δ represents velocity of longitudinal wave situation of change with phase angle in nearly vertical incidence, and p represents that the level is slow
Degree, ε represents the velocity of longitudinal wave in horizontal and vertical difference value.
In one embodiment, the shear wave determining module can include：4th reservoir parameter computing unit, for calculating
The inverse of value of the shear wave velocity in vertical direction in the reservoir to be measured, and square reciprocal is calculated, it is used as the 4th storage
Layer parameter；5th reservoir parameter computing unit, for calculate the velocity of longitudinal wave vertical direction value square with it is described
Product between the biquadratic of horizontal slowness, obtains the 5th reservoir parameter；6th reservoir parameter computing unit, it is described for calculating
Velocity of longitudinal wave vertical direction value and the shear wave velocity in square of the ratio of the value of vertical direction, and calculate described
Ratio square with the horizontal slowness square product as level value, by the ratio square with the level value phase
Multiply, obtain the 6th reservoir parameter；7th reservoir parameter computing unit, for calculating the 5th reservoir parameter and the 6th storage
The difference of layer parameter, is used as the 7th reservoir parameter；8th reservoir parameter computing unit, for calculating the speed difference value and institute
The product of the 7th reservoir parameter is stated, and resulting result of product is added with the 4th reservoir parameter, the knot after will add up
Fruit is used as the 8th reservoir parameter；The vertical slowness computing unit of shear wave, it is slow for calculating the 8th reservoir parameter, the level
Degree square between difference, using resulting difference as the reservoir transverse wave to be measured vertical slowness.
In one embodiment, the shear wave determining module specifically can be used for determining the storage to be measured according to below equation
The vertical slowness of layer plane shear wave：
In above formula, qβRepresent the vertical slowness of the reservoir transverse wave to be measured, vS0Represent the shear wave velocity vertical
The value in direction, vP0The value of velocity of longitudinal wave in the reservoir to be measured in vertical direction is represented, p represents the horizontal slowness, ε tables
Show the velocity of longitudinal wave in horizontal and vertical difference value, δ represents the longitudinal wave velocity degree in nearly vertical incidence with phase angle
Situation of change.
There is provided a kind of vertical slowness formula of the new reservoir plane compressional wave to be measured, institute in embodiments of the present invention
The vertical slowness formula of reservoir transverse wave to be measured is stated, and utilizes the reflectance factor of the VTI media constructed by above-mentioned two formula
The method that calculation carries out AVA invertings to the reservoir to be measured.Above-mentioned two formula not only has sufficiently high computational accuracy,
Also there is relatively simple form and relatively weak non-linear, the reservoir plane compressional wave to be measured proposed using the application
Vertical slowness formula, the vertical slowness formula of the reservoir transverse wave to be measured calculated the reflectance factors of the VTI media
After mode is simplified, when the reflectance factor calculating mode of the VTI media after obtained simplification carries out AVA invertings, good gram
Take in existing inversion technique as it is assumed that the problems such as inversion accuracy is relatively low caused by condition is too many.Further, simplify
The reflectance factor of VTI media afterwards calculates mode and is respectively provided with when strong anisotropy, strong impedance contrast and big incidence angle
Higher computational accuracy.When the inversion result proposed using the application carries out oil field prospecting, surveying accuracy has been obtained effectively
Improve.
Brief description of the drawings
In order to illustrate more clearly of the embodiment of the present application or technical scheme of the prior art,
The required accompanying drawing used in embodiment or description of the prior art will be briefly described below, it is clear that
Ground, drawings in the following description are only some embodiments described in the application, for those of ordinary skill in the art,
Without having to pay creative labor, other accompanying drawings can also be obtained according to these accompanying drawings.
Fig. 1 is a kind of seismic inversion method flow chart that the application is provided；
Fig. 2 is the VTI parameters that a kind of seismic inversion method inverting that the application is provided is obtained, wherein, figure a represents compressional wave speed
Spend vP, scheme b and represent shear wave velocity vS, figure c represents density, schemes the difference results ε that d represents horizontal and vertical p wave interval velocity, and figure e is represented
P wave interval velocity in nearly vertical incidence with phase angle situation of change δ (e)；
Fig. 3 is a kind of structured flowchart for seismic inversion device that the application is provided.
Embodiment
In order that those skilled in the art more fully understand the technical scheme in the application, it is real below in conjunction with the application
The accompanying drawing in example is applied, the technical scheme in the embodiment of the present application is clearly and completely described, it is clear that described implementation
Example only some embodiments of the present application, rather than whole embodiments.Based on the embodiment in the application, this area is common
The every other embodiment that technical staff is obtained under the premise of creative work is not made, should all belong to the application protection
Scope.
When carrying out AVA invertings in view of being based on VTI media in the prior art, due to existing accurate VTI dieletric reflections
The problem of reservoir exploration precision caused by coefficient formula is extremely complex is relatively low, inventors herein proposes a kind of new VTI media anti-
Penetrate the method for simplifying of coefficient calculation, i.e. according to the vertical slowness formula of the reservoir plane P waves to be measured, the storage to be measured
The vertical slowness formula of layer plane SV ripples simplifies to VTI dieletric reflection coefficient formulas, recycles the VTI media after simplifying anti-
Penetrate coefficient calculation and AVA invertings are carried out to reservoir to be measured.Based on this, it is proposed that a kind of seismic inversion method, as shown in figure 1,
It may comprise steps of：
S101：Based on velocity of longitudinal wave in reservoir to be measured in the value of vertical direction, the velocity of longitudinal wave in nearly vertical incidence
When with the situation of change at phase angle, horizontal slowness, the velocity of longitudinal wave in horizontal and vertical difference value, determine the storage to be measured
The vertical slowness relational expression of layer plane compressional wave.
Specifically, the reservoir to be measured can be shale reservoir or sandstone reservoir, carbonate reservoir etc. other
Various types of reservoirs to be measured, the application is not construed as limiting to this.
Wherein, in anisotropic medium, the longitudinal wave propagation speed of different directions is different, and the velocity of longitudinal wave can
To be relevant with incident angle.I.e. it is 0 that the velocity of longitudinal wave in the case of incidence angles degree, which is equal to vertical incidence incidence angle,
Velocity of longitudinal wave when spending is multiplied by an expression formula relevant with angle.
According to the particle vibration direction of ripple and the difference of the direction of propagation, shear wave (S ripples) and compressional wave (P ripples) can be divided into.
Particle vibration direction is referred to as P ripples with direction of propagation identical ripple；The particle vibration direction ripple vertical with the direction of propagation is referred to as S ripples.
Particle vibration occurs ripple in the face perpendicular with the propagation face of ripple and is referred to as shear wave (SV ripples), particle vibration occur with ripple
Ripple in the parallel face in propagation face is referred to as SH ripples.Further, PP ripples can be the back wave after the P ripples incidence, P ripples
Speed can be velocity of longitudinal wave.
In one embodiment of the application, the vertical slowness of the reservoir plane compressional wave to be measured is determined, can be included：Meter
Value of the velocity of longitudinal wave described in the reservoir to be measured in vertical direction is calculated, and calculates square reciprocal of the value, as
First reservoir parameter；Calculate square of the horizontal slowness, and by the horizontal slowness square with the velocity of longitudinal wave near
It is multiplied during vertical incidence with the situation of change at phase angle, is used as the second reservoir parameter；The velocity of longitudinal wave is calculated horizontal and vertical
To difference value, the velocity of longitudinal wave in nearly vertical incidence with phase angle situation of change difference, be used as speed difference value, meter
Square of the velocity of longitudinal wave in the value of vertical direction is calculated, the biquadratic of the horizontal slowness is calculated, the speed difference is calculated
Different value, the value square, the difference between the biquadratic of the horizontal slowness, and resulting difference is used as to the 3rd storage
Layer parameter；The root of the 3rd reservoir parameter is calculated, and regard the root as the vertical of the reservoir plane compressional wave to be measured
Slowness.
Specifically, the vertical slowness of the reservoir plane compressional wave to be measured can be determined according to below equation：
In above formula, qαRepresent the vertical slowness of the reservoir plane compressional wave to be measured, vP0Represent the velocity of longitudinal wave vertical
The value in direction, δ represents velocity of longitudinal wave situation of change with phase angle in nearly vertical incidence, and p represents that the level is slow
Degree, ε represents the velocity of longitudinal wave in horizontal and vertical difference value.
S102：According to shear wave velocity in reservoir to be measured in the value of vertical direction, the velocity of longitudinal wave in vertical direction
Value, the velocity of longitudinal wave are in nearly vertical incidence with the situation of change at phase angle, horizontal slowness, the velocity of longitudinal wave in transverse direction
With vertical difference value, the vertical slowness relational expression of the reservoir transverse wave to be measured is determined.
In one embodiment of the application, the vertical slowness of the reservoir transverse wave to be measured is determined, can be included：Meter
The inverse of value of the shear wave velocity in vertical direction in the reservoir to be measured is calculated, and calculates square reciprocal, the 4th is used as
Reservoir parameter；Calculate square multiplying the biquadratic of the horizontal slowness between of the velocity of longitudinal wave in the value of vertical direction
Product, obtains the 5th reservoir parameter；Value of the velocity of longitudinal wave in vertical direction is calculated with the shear wave velocity in vertical direction
Value ratio square, and calculate the ratio square with the horizontal slowness square product as level value,
By square being multiplied with the level value for the ratio, the 6th reservoir parameter is obtained；Calculate the 5th reservoir parameter and described
The difference of 6th reservoir parameter, is used as the 7th reservoir parameter；Calculate multiplying for the speed difference value and the 7th reservoir parameter
Product, and resulting result of product is added with the 4th reservoir parameter, the result after will add up is used as the 8th reservoir parameter；
Calculate the 8th reservoir parameter, the horizontal slowness square between difference, using resulting difference as described to be measured
The vertical slowness of reservoir transverse wave.
Specifically, determining the vertical slowness of the reservoir transverse wave to be measured according to below equation：
In above formula, qβRepresent the vertical slowness of the reservoir transverse wave to be measured, vS0Represent the shear wave velocity vertical
The value in direction, vP0The value of velocity of longitudinal wave in the reservoir to be measured in vertical direction is represented, p represents the horizontal slowness, ε tables
Show the velocity of longitudinal wave in horizontal and vertical difference value, δ represents the longitudinal wave velocity degree in nearly vertical incidence with phase angle
Situation of change.
S103：According to shear wave in the vertical slowness relational expression of the reservoir plane compressional wave to be measured, the reservoir plane to be measured
Vertical slowness relational expression, build the reflectance factor calculating mode of the new Method in Transverse Isotropic Medium with vertical symmetry axis.
S104：Using the dieletric reflection coefficient calculation, row amplitude is entered with angle change to the reservoir to be measured
Inverting.
In this application, can be based on the Method in Transverse Isotropic Medium with vertical symmetry axis in the reservoir to be measured
Parameter (i.e. VTI parameters) has anisotropic feature, carries out AVA invertings to the reservoir to be measured in the following way：
In this application, the VTI parameters can include：Three elastic parameters and two anisotropic parameterses.Wherein, bullet
Property parameter can be：Velocity of longitudinal wave v of the reservoir to be measured in vertical directionP, reservoir to be measured vertical direction shear wave velocity vSWith it is close
Spend ρ；Anisotropic parameters can be：The difference results ε and p wave interval velocity of horizontal and vertical p wave interval velocity in nearly vertical incidence with
The situation of change δ at phase angle.It is of course also possible to including formation thickness etc. other specification, the application is not construed as limiting to this.
S4-1：Geological data based on the reservoir to be measured extracts angle dependency wavelet, passes through log data forward modeling mould
Fit well lie geological data and determine amplitude scaling factors, and applied to the wavelet；Determine the PP ripples and PS radio frequency channels of angle domain
Collection.
Statistical method is taken to extract wavelet based on actual earthquake prestack trace gather and log data, wavelet is in communication process
Influenceed can occur waveform or frequency change by stratum, amplitude can be effectively improved dependent on the seismic wavelet of incident angle by extracting
With degree.Calculated amplitude zoom factor, and applied to above-mentioned wavelet, reach the amplitude matches of analog record and physical record.Its
In, it is that each seismic channel in trace gather uses unified amplitude scaling factors, to ensure when earthquake data SNR is higher
Amplitude with offset distance variation relation；When signal to noise ratio is relatively low, can divide it is near, in, remote offset distance difference calculated amplitude zoom factor,
Ensure analog record and the best match of physical record, reduce influence of the noise to refutation process.
During due to carrying out inverting individually with the PP ripples in geological data, there is larger uncertainty, and PS ripple packets
Containing abundant density and anisotropic parameters information, therefore, AVA is carried out together by using PP ripples and PS ripples in this application anti-
Drill, so as to improve the stability and precision of inverting.In this application, it can be situated between by combining log data and accurate VTI
Matter reflectance factor equation forward simulation determines the PP ripples and PS radio frequency channel collection of angle domain.
S4-2：Using seismotectonics interpretation data and log data, initial elasticity parameter model is set up.
The model of elastic parameter can be set up using three dimensions interpolation method.First with the method for scatterplot interpolation to each
The data of layer position enter row interpolation, complete geologic horizon modeling, then carry out elastic parameter lateral interpolation according to geologic horizon, will
Well logging information carries out lateral interpolation, so as to calculate the elastic parameter value for obtaining each putting underground, and sets up initial elasticity parameter mould
Type.
S4-3：From the accurate reflection coefficient formula of VTI media, based on weak anisotropy it is assumed that to omit high-order infinite
In a small amount, the reflectance factor for being derived by the new VTI media calculates mode.Based on the initial elasticity parameter model and described
The reflectance factors of VTI media calculates the multi-wave seismic prestack trace gather of mode forward simulation angle domain, and by forward simulation record and
Physical record directly calculates PP ripples and PS ripple inverting residual errors.
By the AVA inversion algorithms that the application is used need to ask for ten unknown elasticity including 5 VTI parameters
The first-order partial derivative of parameter and anisotropic parameters, because accurate VTI reflectance factors solution expression formula is extremely complex, therefore this
Individual first-order partial derivative solves extremely complex.Make discovery from observation, the vertical slowness formula of the reservoir plane P waves to be measured,
The vertical slowness formula of the reservoir PLANE SV WAVES to be measured, the complicated expression-form of the two slownesses directly results in expressed intact
The labyrinth of formula, therefore, present applicant proposes under the premise of computational accuracy loss is less, to the expression formula of the two slownesses
Simplified, making the solution of final first-order partial derivative becomes relatively easy and stably.
In accurate VTI reflectance factors equation, the vertical slowness of different wave modes is represented by：
In above formula,p
=sin (θ)/VP(θ)；
In above formula, qαRepresent the vertical slowness of the reservoir plane P waves to be measured, qβRepresent the reservoir PLANE SV WAVES to be measured
Vertical slowness formula, cijRigidity is represented, ρ represents density, and p represents horizontal slowness, VP(θ) represents P phases velocity of wave, and θ represents phase
Parallactic angle.
Thomsen introduces one group of new anisotropic parameters collection for being applied to VTI media, i.e. can be hung down using two
To speed and three nondimensional anisotropic parameterses rigidity independent to five be replaced：
In above formula, vP0Represent that velocity of longitudinal wave is in the value of vertical direction, v in the reservoir to be measuredS0Represent the storage to be measured
Shear wave velocity is in the value of vertical direction in layer, and γ represents the difference results of horizontal and vertical SH wave velocities, and ε represents p wave interval velocity
In horizontal and vertical difference value, δ describes p wave interval velocity situation of change with phase angle in nearly vertical incidence.Due to
SH ripples are full decoupled from P ripples and SV ripples in VTI media, so only considering the reflection and transmission of P ripples and SV ripples in this application
Coefficient.
Based on weak anisotropy it is assumed that P phases velocity of wave V in VTI mediaP(θ) and rigidity c13Anisotropy can be utilized
Parameter ε and δ linear expression is：
VP(θ)≈vP0(1+δsin2θcos2θ+εsin4θ) (3)
c13≈c33(1+δ)-2c55 (4)
According to equation (2), we can obtain the relatively simple horizontal slowness of form：
Equation (2) and (3) are substituted into the above-mentioned K of equation1、K2、K3Expression formula in, and utilize Thomsen anisotropic parameterses
Independent rigidity is replaced with vertical speed, while omitting higher order term, can be obtained：
By K2And K3Substitute into K1In, it can obtain：
Higher order term dimensionless is omitted, can be obtained：
Order,
Then have
Equation (8) and (9) are substituted into equation (1), can be obtained：
So, the vertical slowness formula of the simple reservoir plane P waves to be measured of form, the storage to be measured can have been obtained
The vertical slowness formula of layer plane SV ripples.
From the accurate reflection coefficient formula of existing VTI media, choose and include hanging down for the plane P waves that the application is proposed
Straight slowness, the accurate reflection coefficient formula of the VTI media of the vertical slowness of PLANE SV WAVES, and described treat form is simple
Vertical slowness formula, the vertical slowness formula of the reservoir PLANE SV WAVES to be measured of reservoir plane P waves, i.e. formula (10) is surveyed to substitute into
In the accurate reflection coefficient formula of the VTI media selected, so as to construct the reflectometer of new VTI media
Calculation mode.
Following formula is the accurate reflection coefficient formula of the VTI media, and the horizontal slowness simplified and vertical slowness are substituted into
Into following formula, it is possible to which realization is approximately simplified to whole equation.
In above formula,
In above formula, RPPRepresent PP wave reflection coefficients, RPSRepresent PS wave reflection coefficients, TPPRepresent PP ripple transmission coefficients, TPSTable
Show PS ripple transmission coefficients, (1) represents top dielectric, and (2) represent layer dielectric.
The reflectance factor based on the initial elasticity parameter model and the VTI media calculates mode again, is answered with reference to above-mentioned
With the wavelet of the amplitude scaling factors, the multi-wave seismic prestack trace gather of forward simulation angle domain.According to above-mentioned forward simulation
Record and physical record, which are directly calculated, obtains PP ripples and PS ripple inverting residual errors.
Afterwards, it is possible to use the reflectance factor of the VTI media after simplifying calculates mode, with reference to shown PP ripples and PS ripples
Inverting residual error carries out AVA invertings, so that it is determined that VTI parameters, i.e.,：Elastic parameter and anisotropic parameters.
S4-4：Log data based on the reservoir to be measured determines the prior density function of VTI parameters.
VTI parameters and its average at some positions of reservoir to be measured is obtained based on the log data, including：Elasticity
Parameter and anisotropic parameters.Assuming that while VTI parameter prior model Gaussian distributeds, also comprising obedience differential La Pula
Vertical piece of bound term of this distribution, obtains the model parameter needed, and calculate five VTI parameters of determination by well logging data analysis
Variance and covariance, build covariance matrix, so as to form the model parameter prior density function R for meeting the work area
(m).It can be expressed as：
In above formula, CmFor the block diagonal matrix comprising elastic parameter and anisotropic parameters correlation, μ is VTI parameters
Mean vector (different VTI parameters need to ask for respectively), D is first order differential operator, klFor scale parameter, l=1,2,3,4,
5, m be the parameter vector to be estimated constituted that formed a line by five parameter vectors.Further, these above-mentioned parameters are not for
Same VTI parameters can choose different values.
S4-5：Object function is rewritten into function on elastic parameter disturbance quantity using generalized linear inversion thought, so
To disturbance quantity derivation and make derivative be equal to the zero iterative formula for obtaining disturbance quantity afterwards, utilize iteration weight weighted least-squares to calculate
Method is solved to disturbance quantity, and the disturbance quantity obtained using solving is updated to model parameter, repeats above job step
Until inverting residual error reaches requirement or reaches maximum iteration, the VTI parameters finally given are exported.
According to Bayesian principles, comprehensive inversion likelihood function and prior density function obtain Posterior probability distribution function.
If VTI parameters are mT=(m1,m2,…,mn)T, observation geological data is dT=(d1,d2,…,dn)T, by Bayes
Theory can be obtained, in the case of known earthquake data before superposition, can be attributed to the problem of inverting underground medium elastic parameter
Solve a posterior probability function：
In above formula, P (d)=∫ p (d | m) p (m) dm are normalization factor, can be regarded as constant, and P (d | m) it is likelihood letter
Number, P (m) is prior probability distribution.
The object function J of inverting is determined according to above-mentioned posterior probability1(m), specifically, can be expressed as follows：
In above formula, dppRepresent actual PP ripples observation geological data, dpsRepresent actual PS ripples observation geological data, Gpp
(m) PP ripple forward modeling earthquake records, G are representedps(m) PS ripple forward modeling earthquake records are represented,Represent that the covariance of PP ripple noises is inverse
Matrix,The covariance inverse matrix of PS ripple noises is represented, R (m) represents prior-constrained item.
Due to the bad direct solution of the model parameter, therefore letter can be carried out to above-mentioned probability function by means of Taylor expansion
Change, become the function on model parameter disturbance quantity, object function is then subjected to derivation to disturbance quantity and makes derivative be equal to zero,
Obtain disturbance quantity Δ m iterative formula.
Δmk=(Hk)-1γk
In above formula,
In above formula,WithRepresent that forward modeling equation is inclined on the single order of ten parameters of interface levels
Derivative Gpp(mk) and Gps(mk) forward modeling earthquake record after kth time iteration is represented respectively.
The VTI dieletric reflections coefficient calculation is substituted into the iterative formula, wherein, the VTI media are anti-
The first-order partial derivative that coefficient accounting equation is penetrated on VTI parameters can be tried to achieve by parsing or numerical method., can be with according to derivation
Obtain the renewal iterative formula of model parameter：
mk+1=mk+ηkΔmk, k=0,1,2 ...
Wherein, ηkIt is the step factor of kth time iteration, m represents the parameter vector being made up of all VTI parameters.In this Shen
Please described in step factor ηk=1, k=4.
PP ripples and PS ripple inverting residual errors, are adopted according to the iterative formula of model parameter disturbance quantity and determined by S4-3
With the disturbance quantity of iteration weight weighted least square algorithm computation model parameter, iteration above-mentioned steps obtain this solution
Disturbance quantity be added on the initial elasticity parameter model, it is possible to obtain new VTI parameters, this new VTI parameter seen
Make new initial model, can be solved again according to the solution expression formula of model parameter disturbance quantity and obtain new disturbance quantity, like this
Repeat to update and iterative, until being finally reached stopping criterion for iteration：PP ripples and PS ripple inverting residual errors determined by S4-3,
And the VTI parameters finally given are exactly final inversion result.
It is illustrated in figure 2 the velocity of longitudinal wave v that a kind of seismic inversion method inverting of the application offer is obtainedP(a), shear wave speed
Spend vS(b), density p (c), the difference results ε (d) of horizontal and vertical p wave interval velocity and p wave interval velocity in nearly vertical incidence with phase
The longitudinal axis represents the time in the situation of change δ (e) at angle, figure, and unit is：Second, transverse axis extremely schemes (e) table respectively from left to right, from figure (a)
Show velocity of longitudinal wave (unit：Km/s), shear wave velocity (unit：Km/s), density (unit：G/cm3), ε and δ.Based on recently like public
The shale reservoir VTI parameter multiwave AVA inversion methods of formula be capable of degree of precision the elastic parameter for predicting VTI media and it is each to
Anisotropic parameter information, the prior distribution of density information is included due to introducing, and carries out joint inversion using PS shear waves, and it is to close
Degree model is also predicted accurately.
Based on same inventive concept, a kind of seismic inversion device is additionally provided in the embodiment of the present invention, such as following implementation
Example is described.Because the principle that seismic inversion device solves problem is similar to seismic inversion method, therefore the reality of seismic inversion device
The implementation that may refer to seismic inversion method is applied, part is repeated and repeats no more.It is used below, term " unit " or " mould
Block " can realize the combination of the software and/or hardware of predetermined function.Although the device described by following examples is preferably with soft
Part is realized, but hardware, or the realization of the combination of software and hardware is also that may and be contemplated.Fig. 3 is of the invention real
A kind of structured flowchart of the seismic inversion device of example is applied, as shown in figure 3, can include：P ripples determining module 301, SV ripples determine mould
Block 302, VTI coefficient determination modules 303, AVA invertings module 304, are illustrated to the structure below.
Compressional wave determining module 301, can be used for the value, described vertical in vertical direction based on velocity of longitudinal wave in reservoir to be measured
Wave velocity is in nearly vertical incidence with the situation of change at phase angle, horizontal slowness, the velocity of longitudinal wave in horizontal and vertical difference
Different value, determines the vertical slowness relational expression of the reservoir plane compressional wave to be measured；
Shear wave determining module 302, can be used for the value, described vertical in vertical direction according to shear wave velocity in reservoir to be measured
Wave velocity when the value of vertical direction, the velocity of longitudinal wave are in nearly vertical incidence with the situation of change at phase angle, horizontal slowness,
The velocity of longitudinal wave determines the vertical slowness relational expression of the reservoir transverse wave to be measured in horizontal and vertical difference value；
Reflectance factor calculates module 303, can be used for vertical slowness formula according to the reservoir plane compressional wave to be measured, institute
The vertical slowness formula of shear wave in reservoir plane to be measured is stated, the new Method in Transverse Isotropic Medium with vertical symmetry axis is built
Reflectance factor calculates mode；
Inverting module 304, can be used for using the reflectance factor calculate mode, the reservoir to be measured is entered row amplitude with
The inverting of angle change.
In one embodiment, the compressional wave determining module can include：First reservoir parameter computing unit, can be used for
Value of the velocity of longitudinal wave described in the reservoir to be measured in vertical direction is calculated, and calculates square reciprocal of the value, is made
For the first reservoir parameter；Second reservoir parameter computing unit, can be used for square for calculating the horizontal slowness, and by the water
Square being multiplied with the velocity of longitudinal wave in nearly vertical incidence with the situation of change at phase angle for flat slowness, is used as the second reservoir ginseng
Number；3rd reservoir parameter computing unit, can be used for calculating the velocity of longitudinal wave in horizontal and vertical difference value, the compressional wave
Speed in nearly vertical incidence with phase angle situation of change difference, as speed difference value, calculate the velocity of longitudinal wave and hanging down
Nogata to value square, calculate the biquadratic of the horizontal slowness, calculate the speed difference value, the value it is flat
Difference between side, the biquadratic of the horizontal slowness, and it regard resulting difference as the 3rd reservoir parameter；Compressional wave it is vertical
Slowness computing unit, can be used for the root for calculating the 3rd reservoir parameter, and regard the root as the reservoir to be measured
The vertical slowness of plane compressional wave.
In one embodiment, the compressional wave determining module specifically can be used for determining the storage to be measured according to below equation
The vertical slowness of layer plane compressional wave：
In above formula, qαRepresent the vertical slowness of the reservoir plane compressional wave to be measured, vP0Represent the velocity of longitudinal wave vertical
The value in direction, δ represents velocity of longitudinal wave situation of change with phase angle in nearly vertical incidence, and p represents that the level is slow
Degree, ε represents the velocity of longitudinal wave in horizontal and vertical difference value.
In one embodiment, the shear wave determining module can include：4th reservoir parameter computing unit, can be used for
The inverse of value of the shear wave velocity in vertical direction in the reservoir to be measured is calculated, and calculates square reciprocal, is used as
Four reservoir parameters；5th reservoir parameter computing unit, can be used for calculating the velocity of longitudinal wave vertical direction value it is flat
The square product between the biquadratic of the horizontal slowness, obtains the 5th reservoir parameter；6th reservoir parameter computing unit, can be with
For calculate the velocity of longitudinal wave vertical direction value and the shear wave velocity the value of vertical direction ratio it is flat
Side, and calculate the ratio square with the horizontal slowness square product as level value, by square of the ratio
It is multiplied with the level value, obtains the 6th reservoir parameter；7th reservoir parameter computing unit, can be used for calculating the 5th storage
The difference of layer parameter and the 6th reservoir parameter, is used as the 7th reservoir parameter；8th reservoir parameter computing unit, can be used for
Calculate the product of the speed difference value and the 7th reservoir parameter, and by resulting result of product and the 4th reservoir
Parameter is added, and the result after will add up is used as the 8th reservoir parameter；The vertical slowness computing unit of shear wave, can be used for calculating institute
State the 8th reservoir parameter, the horizontal slowness square between difference, put down resulting difference as the reservoir to be measured
The vertical slowness of face shear wave.
In one embodiment, the shear wave determining module specifically can be used for determining the storage to be measured according to below equation
The vertical slowness of layer plane shear wave：
In above formula, qβRepresent the vertical slowness of the reservoir transverse wave to be measured, vS0Represent the shear wave velocity vertical
The value in direction, vP0The value of velocity of longitudinal wave in the reservoir to be measured in vertical direction is represented, p represents the horizontal slowness, ε tables
Show the velocity of longitudinal wave in horizontal and vertical difference value, δ represents the longitudinal wave velocity degree in nearly vertical incidence with phase angle
Situation of change.
Although mentioning the vertical slowness determination mode of plane P waves in teachings herein, the vertical slowness of PLANE SV WAVES is determined
Mode and AVA inverting modes etc. are described, and still, the application is not limited to be the feelings described by the embodiment of the present application
Condition.Some professional standards use embodiment amended slightly in self-defined mode or the practice processes of embodiment description
Can also realize above-described embodiment it is identical, equivalent or close or deformation after it is anticipated that implementation result.Using these modifications or change
The vertical slowness determination mode of plane P waves after shape, the vertical slowness determination mode of PLANE SV WAVES and AVA inverting modes etc. are obtained
Within the scope of the embodiment taken, the optional embodiment that still may belong to the application.
Although this application provides the method operating procedure as described in embodiment or flow chart, based on conventional or noninvasive
The means for the property made can include more or less operating procedures.The step of being enumerated in embodiment order is only numerous steps
A kind of mode in execution sequence, unique execution sequence is not represented., can be with when device in practice or end product execution
Performed according to embodiment or method shown in the drawings order or parallel execution (such as parallel processor or multiple threads
Environment, even distributed data processing environment).Term " comprising ", "comprising" or its any other variant are intended to
Nonexcludability is included, so that process, method, product or equipment including a series of key elements not only will including those
Element, but also other key elements including being not expressly set out, or also include being this process, method, product or equipment
Intrinsic key element.In the absence of more restrictions, be not precluded from the process including the key element, method, product or
Also there are other identical or equivalent elements in person's equipment.
Unit, device or module that above-described embodiment is illustrated etc., can specifically be realized by computer chip or entity, or
Realized by the product with certain function.For convenience of description, describe to be divided into various modules point during apparatus above with function
Do not describe.Certainly, when implementing the application can the function of each module in same or multiple softwares and/or hardware it is real
It is existing, the module for realizing same function can also be realized by the combination of multiple submodule or subelement etc..Dress described above
It is only schematical to put embodiment, for example, the division of the unit, only a kind of division of logic function, when actually realizing
There can be other dividing mode, such as multiple units or component can combine or be desirably integrated into another system, or one
A little features can be ignored, or not perform.It is another, shown or discussed coupling or direct-coupling or communication link each other
Connect can be can be electrical, machinery or other shapes by some interfaces, the INDIRECT COUPLING or communication connection of device or unit
Formula.
It is also known in the art that in addition to realizing controller in pure computer readable program code mode, it is complete
Controller can be caused with gate, switch, application specific integrated circuit, programmable by the way that method and step is carried out into programming in logic entirely
Logic controller realizes identical function with the form of embedded microcontroller etc..Therefore this controller is considered one kind
Hardware component, and the device for realizing various functions included to its inside can also be considered as the structure in hardware component.Or
Person even, not only can be able to will be the software module of implementation method but also can be hardware for realizing that the device of various functions is considered as
Structure in part.
The application can be described in the general context of computer executable instructions, such as program
Module.Usually, program module includes performing particular task or realizes routine, program, object, the group of particular abstract data type
Part, data structure, class etc..The application can also be put into practice in a distributed computing environment, in these DCEs,
Task is performed by the remote processing devices connected by communication network.In a distributed computing environment, program module can
With positioned at including in the local and remote computer-readable storage medium including storage device.
As seen through the above description of the embodiments, those skilled in the art can be understood that the application can
Realized by the mode of software plus required general hardware platform.Understood based on such, the technical scheme essence of the application
On the part that is contributed in other words to prior art can be embodied in the form of software product, the computer software product
It can be stored in storage medium, such as ROM/RAM, magnetic disc, CD, including some instructions are to cause a computer equipment
(can be personal computer, mobile terminal, server, or network equipment etc.) performs each embodiment of the application or implementation
Method described in some parts of example.
Each embodiment in this specification is described by the way of progressive, same or analogous portion between each embodiment
Divide mutually referring to what each embodiment was stressed is the difference with other embodiment.The application can be used for crowd
In more general or special purpose computing system environments or configuration.For example：Personal computer, server computer, handheld device or
Portable set, laptop device, multicomputer system, the system based on microprocessor, set top box, programmable electronics are set
Standby, network PC, minicom, DCE of mainframe computer including any of the above system or equipment etc..
Although depicting the application by embodiment, it will be appreciated by the skilled addressee that the application have it is many deformation and
Change is without departing from spirit herein, it is desirable to which appended claim includes these deformations and changed without departing from the application's
Spirit.
Claims (10)
1. a kind of seismic inversion method, it is characterised in that including：
Based on velocity of longitudinal wave in reservoir to be measured when the value of vertical direction, the velocity of longitudinal wave are in nearly vertical incidence with phase angle
Situation of change, horizontal slowness, the velocity of longitudinal wave in horizontal and vertical difference value, determine the reservoir plane compressional wave to be measured
Vertical slowness relational expression；
According to shear wave velocity in reservoir to be measured the value of vertical direction, the velocity of longitudinal wave vertical direction value, described
Velocity of longitudinal wave is in nearly vertical incidence with the situation of change at phase angle, horizontal slowness, the velocity of longitudinal wave horizontal and vertical
Difference value, determines the vertical slowness relational expression of the reservoir transverse wave to be measured；
Closed according to the vertical slowness of the vertical slowness relational expression of the reservoir plane compressional wave to be measured, the reservoir transverse wave to be measured
It is formula, the reflectance factor for building the new Method in Transverse Isotropic Medium with vertical symmetry axis calculates mode；
Mode is calculated using the reflectance factor, inverting of the row amplitude with angle change is entered to the reservoir to be measured.
2. the method as described in claim 1, it is characterised in that determine the vertical slowness of the reservoir plane compressional wave to be measured, bag
Include：
Value of the velocity of longitudinal wave described in the reservoir to be measured in vertical direction is calculated, and calculates the reciprocal of the value and is put down
Side, is used as the first reservoir parameter；
Calculate square of the horizontal slowness, and by the horizontal slowness square with the velocity of longitudinal wave in nearly vertical incidence
It is multiplied with the situation of change at phase angle, is used as the second reservoir parameter；
The velocity of longitudinal wave is calculated when horizontal and vertical difference value, the velocity of longitudinal wave are in nearly vertical incidence with phase angle
The difference of situation of change, as speed difference value, calculates square of value of the velocity of longitudinal wave in vertical direction, calculates the water
The biquadratic of flat slowness, calculate the speed difference value, the value square, the difference between the biquadratic of the horizontal slowness
Value, and it regard resulting difference as the 3rd reservoir parameter；
The root of the 3rd reservoir parameter is calculated, and regard the root as the vertical slow of the reservoir plane compressional wave to be measured
Degree.
3. method as claimed in claim 2, it is characterised in that determine the reservoir plane compressional wave to be measured according to below equation
Vertical slowness：
<mrow>
<msub>
<mi>q</mi>
<mi>&alpha;</mi>
</msub>
<mo>=</mo>
<msqrt>
<mrow>
<mfrac>
<mn>1</mn>
<msubsup>
<mi>v</mi>
<mrow>
<mi>P</mi>
<mn>0</mn>
</mrow>
<mn>2</mn>
</msubsup>
</mfrac>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>+</mo>
<mn>2</mn>
<mi>&delta;</mi>
<mo>)</mo>
</mrow>
<msup>
<mi>p</mi>
<mn>2</mn>
</msup>
<mo>-</mo>
<mn>2</mn>
<mrow>
<mo>(</mo>
<mi>&epsiv;</mi>
<mo>-</mo>
<mi>&delta;</mi>
<mo>)</mo>
</mrow>
<msubsup>
<mi>v</mi>
<mrow>
<mi>P</mi>
<mn>0</mn>
</mrow>
<mn>2</mn>
</msubsup>
<msup>
<mi>p</mi>
<mn>4</mn>
</msup>
</mrow>
</msqrt>
</mrow>
In above formula, qαRepresent the vertical slowness of the reservoir plane compressional wave to be measured, vP0Represent the velocity of longitudinal wave in vertical direction
Value, δ represents velocity of longitudinal wave situation of change with phase angle in nearly vertical incidence, and p represents the horizontal slowness, ε
Represent the velocity of longitudinal wave in horizontal and vertical difference value.
4. the method as described in claim 1, it is characterised in that determine the vertical slowness of the reservoir transverse wave to be measured, bag
Include：
The inverse of value of the shear wave velocity in vertical direction in the reservoir to be measured is calculated, and calculates square reciprocal, is made
For the 4th reservoir parameter；
Product square the biquadratic of the horizontal slowness between of the velocity of longitudinal wave in the value of vertical direction is calculated, is obtained
To the 5th reservoir parameter；
Calculate the velocity of longitudinal wave vertical direction value and the shear wave velocity the value of vertical direction ratio it is flat
Side, and calculate the ratio square with the horizontal slowness square product as level value, by square of the ratio
It is multiplied with the level value, obtains the 6th reservoir parameter；
The difference of the 5th reservoir parameter and the 6th reservoir parameter is calculated, the 7th reservoir parameter is used as；
Calculate the product of the speed difference value and the 7th reservoir parameter, and by resulting result of product and the described 4th
Reservoir parameter is added, and the result after will add up is used as the 8th reservoir parameter；
Calculate the 8th reservoir parameter, the horizontal slowness square between difference, using resulting difference as described
The vertical slowness of reservoir transverse wave to be measured.
5. method as claimed in claim 4, it is characterised in that determine the reservoir transverse wave to be measured according to below equation
Vertical slowness：
<mrow>
<msub>
<mi>q</mi>
<mi>&beta;</mi>
</msub>
<mo>=</mo>
<msqrt>
<mrow>
<mfrac>
<mn>1</mn>
<msubsup>
<mi>v</mi>
<mrow>
<mi>S</mi>
<mn>0</mn>
</mrow>
<mn>2</mn>
</msubsup>
</mfrac>
<mo>+</mo>
<mn>2</mn>
<mrow>
<mo>(</mo>
<mi>&epsiv;</mi>
<mo>-</mo>
<mi>&delta;</mi>
<mo>)</mo>
</mrow>
<mrow>
<mo>(</mo>
<msubsup>
<mi>v</mi>
<mrow>
<mi>P</mi>
<mn>0</mn>
</mrow>
<mn>2</mn>
</msubsup>
<msup>
<mi>p</mi>
<mn>4</mn>
</msup>
<mo>-</mo>
<mfrac>
<msubsup>
<mi>v</mi>
<mrow>
<mi>P</mi>
<mn>0</mn>
</mrow>
<mn>2</mn>
</msubsup>
<msubsup>
<mi>v</mi>
<mrow>
<mi>S</mi>
<mn>0</mn>
</mrow>
<mn>2</mn>
</msubsup>
</mfrac>
<msup>
<mi>p</mi>
<mn>2</mn>
</msup>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msup>
<mi>p</mi>
<mn>2</mn>
</msup>
</mrow>
</msqrt>
</mrow>
In above formula, qβRepresent the vertical slowness of the reservoir transverse wave to be measured, vS0Represent the shear wave velocity in vertical direction
Value, vP0The value of velocity of longitudinal wave in the reservoir to be measured in vertical direction is represented, p represents the horizontal slowness, and ε represents institute
Velocity of longitudinal wave is stated in horizontal and vertical difference value, δ represents longitudinal wave velocity degree change with phase angle in nearly vertical incidence
Change situation.
6. a kind of seismic inversion device, it is characterised in that including：
Compressional wave determining module, for based on velocity of longitudinal wave in reservoir to be measured in the value of vertical direction, the velocity of longitudinal wave near
During vertical incidence institute is determined with the situation of change at phase angle, horizontal slowness, the velocity of longitudinal wave in horizontal and vertical difference value
State the vertical slowness relational expression of reservoir plane compressional wave to be measured；
Shear wave determining module, for being hung down according to shear wave velocity in reservoir to be measured in the value of vertical direction, the velocity of longitudinal wave
Nogata to value, the velocity of longitudinal wave in nearly vertical incidence with the situation of change at phase angle, horizontal slowness, the compressional wave speed
Degree determines the vertical slowness relational expression of the reservoir transverse wave to be measured in horizontal and vertical difference value；
Reflectance factor calculates module, for the vertical slowness relational expression according to the reservoir plane compressional wave to be measured, the storage to be measured
The vertical slowness relational expression of shear wave in layer plane, builds the reflection system of the new Method in Transverse Isotropic Medium with vertical symmetry axis
Number calculation；
Inverting module, for calculating mode using the reflectance factor, row amplitude is entered with angle change to the reservoir to be measured
Inverting.
7. device as claimed in claim 6, it is characterised in that the compressional wave determining module includes：
First reservoir parameter computing unit, for calculating value of the velocity of longitudinal wave described in the reservoir to be measured in vertical direction,
And square reciprocal of the value is calculated, it is used as the first reservoir parameter；
Second reservoir parameter computing unit, square for calculating the horizontal slowness, and by the horizontal slowness square with
The velocity of longitudinal wave is multiplied in nearly vertical incidence with the situation of change at phase angle, is used as the second reservoir parameter；
3rd reservoir parameter computing unit, for calculating the velocity of longitudinal wave in horizontal and vertical difference value, compressional wave speed
The difference of the situation of change in nearly vertical incidence with phase angle is spent, as speed difference value, the velocity of longitudinal wave is calculated vertical
Square of the value in direction, calculates the biquadratic of the horizontal slowness, calculate the speed difference value, the value square,
Difference between the biquadratic of the horizontal slowness, and it regard resulting difference as the 3rd reservoir parameter；
The vertical slowness computing unit of compressional wave, the root for calculating the 3rd reservoir parameter, and it regard the root as institute
State the vertical slowness of reservoir plane compressional wave to be measured.
8. device as claimed in claim 7, it is characterised in that the compressional wave determining module is specifically for true according to below equation
The vertical slowness of the fixed reservoir plane compressional wave to be measured：
<mrow>
<msub>
<mi>q</mi>
<mi>&alpha;</mi>
</msub>
<mo>=</mo>
<msqrt>
<mrow>
<mfrac>
<mn>1</mn>
<msubsup>
<mi>v</mi>
<mrow>
<mi>P</mi>
<mn>0</mn>
</mrow>
<mn>2</mn>
</msubsup>
</mfrac>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>+</mo>
<mn>2</mn>
<mi>&delta;</mi>
<mo>)</mo>
</mrow>
<msup>
<mi>p</mi>
<mn>2</mn>
</msup>
<mo>-</mo>
<mn>2</mn>
<mrow>
<mo>(</mo>
<mi>&epsiv;</mi>
<mo>-</mo>
<mi>&delta;</mi>
<mo>)</mo>
</mrow>
<msubsup>
<mi>v</mi>
<mrow>
<mi>P</mi>
<mn>0</mn>
</mrow>
<mn>2</mn>
</msubsup>
<msup>
<mi>p</mi>
<mn>4</mn>
</msup>
</mrow>
</msqrt>
</mrow>
In above formula, qαRepresent the vertical slowness of the reservoir plane compressional wave to be measured, vP0Represent the velocity of longitudinal wave in vertical direction
Value, δ represents velocity of longitudinal wave situation of change with phase angle in nearly vertical incidence, and p represents the horizontal slowness, ε
Represent the velocity of longitudinal wave in horizontal and vertical difference value.
9. device as claimed in claim 6, it is characterised in that the shear wave determining module includes：
4th reservoir parameter computing unit, falls for calculating shear wave velocity in the reservoir to be measured in the value of vertical direction
Number, and square reciprocal is calculated, it is used as the 4th reservoir parameter；
5th reservoir parameter computing unit, for calculate the velocity of longitudinal wave vertical direction value square with the level
Product between the biquadratic of slowness, obtains the 5th reservoir parameter；
6th reservoir parameter computing unit, exists for calculating the velocity of longitudinal wave in the value of vertical direction with the shear wave velocity
Square of the ratio of the value of vertical direction, and calculate the ratio square with the horizontal slowness square product conduct
Level value, by square being multiplied with the level value for the ratio, obtains the 6th reservoir parameter；
7th reservoir parameter computing unit, the difference for calculating the 5th reservoir parameter and the 6th reservoir parameter is made
For the 7th reservoir parameter；
8th reservoir parameter computing unit, the product for calculating the speed difference value and the 7th reservoir parameter, and will
Resulting result of product is added with the 4th reservoir parameter, and the result after will add up is used as the 8th reservoir parameter；
The vertical slowness computing unit of shear wave, for calculate the 8th reservoir parameter, the horizontal slowness square between
Difference, using resulting difference as the reservoir transverse wave to be measured vertical slowness.
10. device as claimed in claim 9, it is characterised in that the shear wave determining module is specifically for according to below equation
Determine the vertical slowness of the reservoir transverse wave to be measured：
<mrow>
<msub>
<mi>q</mi>
<mi>&beta;</mi>
</msub>
<mo>=</mo>
<msqrt>
<mrow>
<mfrac>
<mn>1</mn>
<msubsup>
<mi>v</mi>
<mrow>
<mi>S</mi>
<mn>0</mn>
</mrow>
<mn>2</mn>
</msubsup>
</mfrac>
<mo>+</mo>
<mn>2</mn>
<mrow>
<mo>(</mo>
<mi>&epsiv;</mi>
<mo>-</mo>
<mi>&delta;</mi>
<mo>)</mo>
</mrow>
<mrow>
<mo>(</mo>
<msubsup>
<mi>v</mi>
<mrow>
<mi>P</mi>
<mn>0</mn>
</mrow>
<mn>2</mn>
</msubsup>
<msup>
<mi>p</mi>
<mn>4</mn>
</msup>
<mo>-</mo>
<mfrac>
<msubsup>
<mi>v</mi>
<mrow>
<mi>P</mi>
<mn>0</mn>
</mrow>
<mn>2</mn>
</msubsup>
<msubsup>
<mi>v</mi>
<mrow>
<mi>S</mi>
<mn>0</mn>
</mrow>
<mn>2</mn>
</msubsup>
</mfrac>
<msup>
<mi>p</mi>
<mn>2</mn>
</msup>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msup>
<mi>p</mi>
<mn>2</mn>
</msup>
</mrow>
</msqrt>
</mrow>
In above formula, qβRepresent the vertical slowness of the reservoir transverse wave to be measured, vS0Represent the shear wave velocity in vertical direction
Value, vP0The value of velocity of longitudinal wave in the reservoir to be measured in vertical direction is represented, p represents the horizontal slowness, and ε represents institute
Velocity of longitudinal wave is stated in horizontal and vertical difference value, δ represents longitudinal wave velocity degree change with phase angle in nearly vertical incidence
Change situation.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710280917.5A CN107092029B (en) | 2017-04-26 | 2017-04-26 | A kind of seismic inversion method and device |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710280917.5A CN107092029B (en) | 2017-04-26 | 2017-04-26 | A kind of seismic inversion method and device |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107092029A true CN107092029A (en) | 2017-08-25 |
| CN107092029B CN107092029B (en) | 2019-01-18 |
Family
ID=59638002
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710280917.5A Active CN107092029B (en) | 2017-04-26 | 2017-04-26 | A kind of seismic inversion method and device |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107092029B (en) |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110187384A (en) * | 2019-06-19 | 2019-08-30 | 湖南科技大学 | Bayes's time-lapse seismic difference inversion method and device |
| CN110261896A (en) * | 2019-04-26 | 2019-09-20 | 中国石油化工股份有限公司 | Stable anisotropy TI medium the Forward Modeling |
| CN110703318A (en) * | 2019-09-24 | 2020-01-17 | 自然资源部第一海洋研究所 | Direct inversion method of pre-stack seismic data |
| CN111208560A (en) * | 2020-01-15 | 2020-05-29 | 中南大学 | Synchronous prediction method for horizontal fractures and vertical fractures of orthogonal medium fractured reservoir |
| CN111239805A (en) * | 2020-02-13 | 2020-06-05 | 中国石油大学(北京) | Block constraint time-lapse seismic difference inversion method and system based on reflectivity method |
| CN112363224A (en) * | 2020-10-29 | 2021-02-12 | 中国石油天然气集团有限公司 | Method and device for judging application potential of transverse wave seismic exploration |
| CN116148925A (en) * | 2023-04-18 | 2023-05-23 | 山东省科学院海洋仪器仪表研究所 | VTI medium spherical longitudinal wave reflection coefficient analysis method |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102012524A (en) * | 2010-09-29 | 2011-04-13 | 中国海洋石油总公司 | Quantitative assessment method for feathering drifting of maritime three-dimensional seismological observation system |
| CN103827696A (en) * | 2011-07-27 | 2014-05-28 | 普拉德研究及开发股份有限公司 | Multi-well anisotropy inversion |
| CN104570085A (en) * | 2013-10-29 | 2015-04-29 | 中国石油化工股份有限公司 | Longitudinal and transverse wave ray parameter domain joint inversion method |
| CN106556861A (en) * | 2015-09-24 | 2017-04-05 | 中国石油化工股份有限公司 | A kind of azimuthal AVO inversion method based on Omnibearing earthquake auto data |
-
2017
- 2017-04-26 CN CN201710280917.5A patent/CN107092029B/en active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102012524A (en) * | 2010-09-29 | 2011-04-13 | 中国海洋石油总公司 | Quantitative assessment method for feathering drifting of maritime three-dimensional seismological observation system |
| CN103827696A (en) * | 2011-07-27 | 2014-05-28 | 普拉德研究及开发股份有限公司 | Multi-well anisotropy inversion |
| CN104570085A (en) * | 2013-10-29 | 2015-04-29 | 中国石油化工股份有限公司 | Longitudinal and transverse wave ray parameter domain joint inversion method |
| CN106556861A (en) * | 2015-09-24 | 2017-04-05 | 中国石油化工股份有限公司 | A kind of azimuthal AVO inversion method based on Omnibearing earthquake auto data |
Non-Patent Citations (4)
| Title |
|---|
| LEON THOMSEN 等: "Weak elastic anisotropy in global seismology", 《THE GEOLOGICAL SOCIETY OF AMERICA SPECIAL PAPER》 * |
| O.PEDERSEN 等: "Wide-angle phase-slowness approximations in VTI media", 《EAGE 68TH CONFERENCE & EXHIBITION》 * |
| 何现启: "EDA介质中地震波传播特征及参数反演研究", 《中国博士学位论文全文数据库 基础科学辑》 * |
| 李景叶 等: "盖层各向异性对时移地震AVO的影响", 《石油物探》 * |
Cited By (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110261896A (en) * | 2019-04-26 | 2019-09-20 | 中国石油化工股份有限公司 | Stable anisotropy TI medium the Forward Modeling |
| CN110261896B (en) * | 2019-04-26 | 2021-07-20 | 中国石油化工股份有限公司 | Stable anisotropic TI medium forward modeling method |
| CN110187384A (en) * | 2019-06-19 | 2019-08-30 | 湖南科技大学 | Bayes's time-lapse seismic difference inversion method and device |
| CN110187384B (en) * | 2019-06-19 | 2021-07-20 | 湖南科技大学 | Bayes time-lapse seismic difference inversion method and device |
| CN110703318A (en) * | 2019-09-24 | 2020-01-17 | 自然资源部第一海洋研究所 | Direct inversion method of pre-stack seismic data |
| CN110703318B (en) * | 2019-09-24 | 2021-06-11 | 自然资源部第一海洋研究所 | Direct inversion method of pre-stack seismic data |
| CN111208560A (en) * | 2020-01-15 | 2020-05-29 | 中南大学 | Synchronous prediction method for horizontal fractures and vertical fractures of orthogonal medium fractured reservoir |
| CN111208560B (en) * | 2020-01-15 | 2021-04-23 | 中南大学 | Synchronous prediction method for horizontal fractures and vertical fractures of orthogonal medium fractured reservoir |
| CN111239805A (en) * | 2020-02-13 | 2020-06-05 | 中国石油大学(北京) | Block constraint time-lapse seismic difference inversion method and system based on reflectivity method |
| CN111239805B (en) * | 2020-02-13 | 2021-02-05 | 中国石油大学(北京) | Block constraint time-lapse seismic difference inversion method and system based on reflectivity method |
| CN112363224A (en) * | 2020-10-29 | 2021-02-12 | 中国石油天然气集团有限公司 | Method and device for judging application potential of transverse wave seismic exploration |
| CN116148925A (en) * | 2023-04-18 | 2023-05-23 | 山东省科学院海洋仪器仪表研究所 | VTI medium spherical longitudinal wave reflection coefficient analysis method |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107092029B (en) | 2019-01-18 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107092029A (en) | A kind of seismic inversion method and device | |
| US10816686B2 (en) | Seismic constrained discrete fracture network | |
| Zhang et al. | Parameter prediction of hydraulic fracture for tight reservoir based on micro-seismic and history matching | |
| CN104375188B (en) | Seismic wave transmission attenuation compensation method and device | |
| US10324211B2 (en) | Seismic spectral balancing | |
| US10359529B2 (en) | Singularity spectrum analysis of microseismic data | |
| US9659252B2 (en) | Method to characterize heterogeneous anisotropic media | |
| US10670754B2 (en) | System and method for processing microseismic data | |
| Cho et al. | Generalized multiscale finite elements for simulation of elastic-wave propagation in fractured media | |
| CN109655894B (en) | Construction method and system of carbonate rock ancient river channel seismic inversion low-frequency model | |
| AU2008239658A1 (en) | Inverse-vector method for smoothing dips and azimuths | |
| CN110501744A (en) | Hydrocarbon source rock organic carbon geophysics quantitative forecasting technique, device, equipment and storage medium | |
| Hong et al. | On a wavelet-based method for the numerical simulation of wave propagation | |
| Genovese et al. | Energy‐compatible modulating functions for the stochastic generation of fully non‐stationary artificial accelerograms and their effects on seismic site response analysis | |
| CN107179545A (en) | The method and apparatus of Nonlinear A VO invertings | |
| CN112147677A (en) | Oil and gas reservoir parameter tag data generation method and device | |
| GB2523460A (en) | Singularity spectrum analysis of microseismic data | |
| CN111025393B (en) | Reservoir prediction method, device, equipment and medium for stratum containing thin coal seam | |
| CN107976711A (en) | A kind of earthquake information quantitative forecast hydrocarbon source rock organic carbon content method and device | |
| Bayuk | Why anisotropy is important for location of microearthquake events in shale? | |
| CN104062680B (en) | A kind of method calculating wave impedance inversion target function gradient | |
| CN114706126B (en) | Quantitative carving method and system for seismic geologic body under geological awareness constraint | |
| CN104297781A (en) | Ray parameter domain unconstrained elastic parameter inversion method | |
| US20240210580A1 (en) | Near Surface Modeling and Drilling Hazard Identification for a Subterranean Formation | |
| Zheng et al. | Improved numerical solution of anisotropic poroelastic wave equation in microseismicity: Graphic process unit acceleration and moment tensor implementation |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |