CN108828660B - A kind of fracture medium PP wave and division PS wave AVO joint inversion method - Google Patents
A kind of fracture medium PP wave and division PS wave AVO joint inversion method Download PDFInfo
- Publication number
- CN108828660B CN108828660B CN201810894632.5A CN201810894632A CN108828660B CN 108828660 B CN108828660 B CN 108828660B CN 201810894632 A CN201810894632 A CN 201810894632A CN 108828660 B CN108828660 B CN 108828660B
- Authority
- CN
- China
- Prior art keywords
- wave
- reflection coefficient
- model
- vertical speed
- inverting
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 43
- 239000011159 matrix material Substances 0.000 claims description 28
- 238000006243 chemical reaction Methods 0.000 claims description 6
- 239000002131 composite material Substances 0.000 claims description 6
- 239000011435 rock Substances 0.000 claims description 5
- 238000004458 analytical method Methods 0.000 claims description 4
- 238000011161 development Methods 0.000 claims description 3
- 238000001615 p wave Methods 0.000 claims description 3
- 230000005570 vertical transmission Effects 0.000 claims description 3
- 238000007906 compression Methods 0.000 abstract description 6
- 230000006835 compression Effects 0.000 abstract description 6
- 230000000694 effects Effects 0.000 abstract description 6
- 239000012530 fluid Substances 0.000 abstract description 3
- 238000012360 testing method Methods 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 2
- 238000010586 diagram 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
- 241000208340 Araliaceae Species 0.000 description 1
- 239000004215 Carbon black (E152) Substances 0.000 description 1
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 1
- 235000003140 Panax quinquefolius Nutrition 0.000 description 1
- 230000002411 adverse Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 230000021615 conjugation Effects 0.000 description 1
- 238000013144 data compression Methods 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 235000008434 ginseng Nutrition 0.000 description 1
- 229930195733 hydrocarbon Natural products 0.000 description 1
- 150000002430 hydrocarbons Chemical class 0.000 description 1
- 230000035945 sensitivity 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
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 fracture medium PP waves and division PS wave AVO joint inversion method, prestack AVA angle gathers have carried out time compression fit, reflection coefficient approximate formula based on R ü ger about HTI medium again, with Gaussian weighting marks algorithm, carry out PP wave, PS1 wave and PS2 wave joint inversion, to obtain the anisotropic parameters and elastic parameter result of fracture medium, suitable for crack HTI medium and isotropic medium, this is because in reflection coefficient approximate formula, when anisotropic parameters are 0, as isotropic medium, even if initial model is linear, inversion method of the invention can also obtain good effect.The multiwave AVO inversion method of HTI medium can provide the effective informations such as lithology and fluid and fracture spacing for fracture medium stratum.
Description
Technical field
The present invention relates to seismic exploration technique fields, and in particular to a kind of fracture medium PP wave is combined with division PS wave AVO
Inversion method.
Background technique
AVO(amplitude variation with offset)/AVA(amplitude variation with
Angle) inversion technique is based on elastic wave theory, using prestack seismic gather reflected amplitude with geophone offset (or incidence angle)
Changing rule, to estimate the physical parameter of predicting reservoir and the important technology of rock characteristic.When seismic wave is incident on reflection circle
When face, Zoeppritz equation describes the pass between the incidence of longitudinal wave incidence shear wave and the reflection and transmission longitudinal wave shear wave of each self-excitation
System, provides theoretical basis for AVO technology；Aki and Richards simplifies Zoeppritz equation, has obtained based on longitudinal and shear wave
The PP wave and PS wave AVO approximate formula of speed and density；Goodway analyzes Lame constants to the sensitivity of hydrocarbon
AVO approximate formula is derived；Gray derived on the formula based on Aki and Richards expressed with Lame constants and density it is anti-
Penetrate coefficient approximate equation；R ü ger has derived the AVO reflection approximate formula expressed with Thomsen anisotropic parameters, and to its public affairs
The accuracy and Applicable media of formula are compared and have been analyzed.
The HTI medium of crack-induced is that the crack of one group of vertical arrangement is embedded in one of isotropic medium theory
EFFECTIVE MEDIUM, the medium are the anisotropic medium with transversely and horizontally symmetry axis.In HTI medium, AVO inverting is carried out to it
More parameters are needed, this is because the anisotropy of its medium can cause reflection coefficient to change, AVO inverting also needs to consider to be situated between
The anisotropic parameters of matter.AVO inverting research for anisotropic medium, due to the more and complexity of its parameter, with previous
AVO inversion method tend not to solve its multi-solution, and then good prediction effect cannot be obtained.The it is proposeds such as Lee part
The method of minimal solution is easier to solve the problems, such as anisotropic medium；Auger etc., Verie and Landro all propose PP wave and
The joint inversion of PS wave solves more Xie Wenti, and obtains preferable achievement in isotropic medium；Fatti etc. utilizes weighting
The accurately inverting P-wave And S impedance of the method for superposition；Nadri and Hartly is with conjugation iterative algorithm inverting elastic parameter and respectively
Anisotropy parameter；Lu etc. is to VTI medium more wave joint inversions anisotropic parameters and elastic parameter 5.
Summary of the invention
In view of the deficiencies of the prior art, the present invention is intended to provide a kind of fracture medium PP wave is combined instead with division PS wave AVO
Method is drilled, time compression fit has been carried out to prestack AVA angle gathers, then approximate about the reflection coefficient of HTI medium based on R ü ger
Formula carries out PP wave, PS1 wave and PS2 wave joint inversion, to obtain each of fracture medium with Gaussian weighting marks algorithm
Anisotropy parameter and elastic parameter result.
To achieve the goals above, the present invention adopts the following technical scheme:
A kind of fracture medium PP wave and division PS wave AVO joint inversion method, include the following steps:
S1, PP wave AVA trace gather, PS1 wave AVA trace gather and PS2 wave AVA trace gather are extracted, and the time of PS1 wave, PS2 wave is pressed
It is reduced on the time of PP wave；
The foundation of S2, initial model: the vertical speed of initial longitudinal wave, the vertical speed of fast transverse wave are obtained by log data
Degree and density model, and according to rock physics well logging information or velocity analysis, establish initial anisotropic parameters model；
S3, inverting:
S3.1 carries out the vertical speed β and density p of isotropic elasticity parameter, that is, longitudinal wave vertical speed α, fast transverse wave
Inverting:
S3.1.1, according to the pattern number of the vertical speed of initial longitudinal wave, the vertical speed of fast transverse wave and density model
According to, the reflection coefficient of PP wave, PS1 wave, PS2 wave is calculated separately with the reflection coefficient approximate formula of R ü ger, at this time anisotropy join
Number ε(V)、δ(V), γ 0, not to anisotropic parameters carry out inverting；The reflection coefficient approximate formula of the R ü ger is as follows:
Wherein, Rpp、Rps1、Rps2The reflection coefficient of PP wave, PS1 wave, PS2 wave is respectively indicated, α indicates the vertical speed of longitudinal wave
Degree, β indicate the vertical speed of fast transverse wave, ε(V)Refer to the difference of longitudinal wave horizontal velocity and vertical speed and the ratio of longitudinal wave vertical speed
Value, indicates the anisotropic degree of longitudinal wave；δ(V)Indicate longitudinal wave and shear wave anisotropy relative size；γ indicates shear wave velocity
Anisotropy development degree；Z indicates that vertical p-wave impedance Z=ρ α, G indicate shear modulus G=ρ of the wave of corresponding vertical transmission
β2, i is incident compressional angle, and j is fast transverse wave angle of reflection, and θ is incident orientation angle, and subscript 1 represents the physical parameter of top dielectric, under
Mark 2 represents the parameter of layer dielectric；
S3.1.2, using theoretical Ricker wavelet, the reflection coefficient difference of calculated PP wave, PS1 wave, PS2 wave
Make the composite traces of PP wave, PS1 wave, PS2 wave；
S3.1.3, joint PP wave, PS1 wave, PS2 wave calculating target function are as follows:
Q=| | (H+ λ I) Δ Μ-JT(Robs-RM)||2；
Wherein, M=(α β ρ ε δ γ)TFor parameter model matrix, according to gauss-newton method, near true value R into
Row first order Taylor is unfolded to obtain J Δ M=Rk+1-Rk；K represents the number of iterations, Δ M=(Δ α Δ β Δ ρ Δ ε Δ δ Δ γ)T
It is the model updated, J is the Jacobian matrix of HTI medium；RobsAnd RMRespectively observe reflection coefficient matrix and model reflection system
Matrix number；
S3.1.4, when objective function is unsatisfactory for required precision, then modify the vertical speed of longitudinal wave, the vertical speed of fast transverse wave
And the model data of density model, and the reflecting system of return step S3.1.1 calculates step again, carries out inverting again；When
Objective function meets required precision, stops iteration；Finally obtain the vertical speed for reaching the inversion result longitudinal wave of well logging resolution ratio
α, the vertical speed β of fast transverse wave and density p；
S3.2, using the inversion result of step S3.1 as constraint condition, to anisotropic parameters carry out inverting:
S3.2.1, using the model data of initial anisotropic parameters model, with the reflection coefficient approximate formula of R ü ger
Calculate separately the reflection coefficient of PP wave, PS1 wave, PS2 wave；
S3.2.2, using theoretical Ricker wavelet, the reflection coefficient difference of calculated PP wave, PS1 wave, PS2 wave
Make the composite traces of PP wave, PS1 wave, PS2 wave；
S3.2.3, combine PP wave, PS1 wave, PS2 wave calculating target function according to the method for step S3.1.3；
S3.2.4, when objective function is unsatisfactory for required precision, then modify the vertical speed of longitudinal wave, the vertical speed of fast transverse wave
And the model data of density model, and the reflecting system of return step S3.2.1 calculates step again, carries out inverting again；When
Objective function meets required precision, stops iteration；Finally obtain the inversion result anisotropic parameters for reaching well logging resolution ratio
ε(V)、δ(V)、γ。
Further, the detailed process of step S1 are as follows:
S1.1, PP wave AVA trace gather extract；
1.1.1 the depth of reflection point, and the corresponding incidence angle of different shot points) are found out according to the trace header information of seismic channel；
1.1.2) for needing the incidence angle of inverting, its amplitude value is extracted, new PP wave AVA trace gather is formed；
S1.2, PS1 wave AVA trace gather extract:
1.2.1 the depth of reflection point, and the corresponding incidence angle of different shot points) are found out according to the trace header information of seismic channel；
1.2.2) for needing the incidence angle of inverting, its amplitude value is extracted, new PS1 wave AVA trace gather is formed；
1.2.3) time of PS1 wave is compressed on the vertical hourage of PP wave by joint velocity of longitudinal wave, after obtaining compression
PS1 wave AVA trace gather:
tPS1(tPP)=tPP+u(tPP)；
tppIndicate the vertical hourage of PP wave, tps1(tpp) indicate PS1 wave time；u(tPP) it is the conversion time difference；
S1.3, PS2 wave AVA trace gather extract
S1.3.1, it is based onCalculate the speed of slow shear-wave；Wherein, γ indicates the anisotropy hair of shear wave velocity
Educate degree, Vs2Indicate the speed of slow shear-wave, Vs1Indicate the speed of fast transverse wave；It is asked further according to the information of the PS1 wave in step S1.2
The depth of reflection point out, and the corresponding incidence angle of different shot points；
S1.3.2, the incidence angle for needing inverting, extract its amplitude value, form new PS2 wave AVA trace gather；
The time of PS2 wave is compressed on the vertical hourage of PP wave by S1.3.3, joint velocity of longitudinal wave, after obtaining compression
PS2 wave AVA trace gather:
tps2(tpp)=tpp+u(tpp)；
tppIndicate the vertical hourage of PP wave, tps2(tpp) indicate PS2 wave time；u(tPP) it is the conversion time difference.
Further, in the reflection coefficient approximate formula of the R ü ger:
Further, in the objective function, extension expression formula of the Jacobian matrix J in AVO are as follows:
Single order local derviation of the reflection coefficient about elastic parameter and anisotropic parameters, is indicated with centered Finite Difference Methods:
With PP wave, PS1 wave, PS2 wave AVA joint inversion, then Jacobian matrix J are as follows:
And H=JTJ then has:
Wherein, ePP、ePS1And ePS2It is PP wave respectively, the residual error of PS1 wave, PS2 wave AVO is all noises and other
The summation of error.
The beneficial effects of the present invention are:
The present invention has carried out time compression fit, then the reflection based on R ü ger about HTI medium to prestack AVA angle gathers
Coefficient approximate formula carries out PP wave, PS1 wave and PS2 wave joint inversion, to obtain crack with Gaussian weighting marks algorithm
Anisotropic parameters and the elastic parameter of medium are as a result, be suitable for crack HTI medium and isotropic medium, this is because reflection
In coefficient approximate formula, when anisotropic parameters are 0, as isotropic medium, even if initial model is linear, the present invention
Inversion method can also obtain good effect.The multiwave AVO inversion method of HTI medium can provide for fracture medium stratum
The effective informations such as lithology and fluid and fracture spacing.
Detailed description of the invention
Fig. 1 is the AVA trace gather schematic diagram of PP wave, PS1 wave and PS2 wave in the test of the embodiment of the present invention；
Fig. 2 is inversion result schematic diagram in the test of the embodiment of the present invention.
Specific embodiment
Below with reference to attached drawing, the invention will be further described, it should be noted that the present embodiment is with this technology side
The premise of case, the detailed implementation method and specific operation process are given, but protection scope of the present invention is not limited to this reality
Apply example.
First the principle of method provided in this embodiment is illustrated and described below.
1, reflection coefficient
Since in the HTI medium of crack, separating phenomenon can occur for shear wave, fast parallel with fracture surface of direction of vibration is split into
Shear wave and the direction of vibration slow shear-wave vertical with fracture surface, thus to crack HTI medium, with PP, PS1, PS2 wave reflection system
Several display expression formulas is conducive to the inverting of AVO.In the present embodiment, HTI medium is carried out using the approximate formula of R ü ger
Inverting:
Wherein, α indicates the vertical speed of longitudinal wave, and β indicates the vertical speed of fast transverse wave, ε(V)Refer to longitudinal wave horizontal velocity and hangs down
The difference of straight speed and the ratio of longitudinal wave vertical speed, indicate the anisotropic degree of longitudinal wave；δ(V)Indicate longitudinal wave and shear wave respectively to
Anisotropic relative size；The anisotropy development degree of γ expression shear wave velocity.Wherein Z indicates vertical p-wave impedance Z=ρ α, G table
Show shear modulus G=ρ β of the wave of corresponding vertical transmission2, i is incident compressional angle, and j is fast transverse wave angle of reflection, and θ is incidence side
Parallactic angle, it is longitudinal wave that subscript PP, PS1, PS2, which respectively refer to incidence wave, and back wave is the plane wave type of longitudinal wave, fast transverse wave, slow shear-wave.
Subscript 1 represents the physical parameter of top dielectric, and subscript 2 represents the parameter of layer dielectric, in above formula also:
According to reflection coefficient, have with the when window of n sampled point come inverting, reflection coefficient matrix expression:
Here θ1…θiRepresent each incidence angle of longitudinal wave, R1…RiRepresent PP wave corresponding to each angle,
The reflection coefficient of PS1 wave, PS2 wave.PP wave, PS1 wave, PS2 wave AVA seismic data as a result, are as follows:
S=WR, (6)
W is Wavelet Martrix, description are as follows:
For HTI medium A VA inverting, nonlinear least square problem is solved with Gauss-Newton method.Given observation
Data B=(b1…bm) and n group model parameter X=(x1…xn), by initial guess model data X(0)Join as first time iteration
Number, Gauss-Newton Methods can make model parameter approaching to reality value by iterating to calculate minimum variance Q.
2, model iteration
What is used in R ü ger (2002) approximate formula is the model parameter of elastic contrastHowever this implementation
Want to go out each layer attribute, therefore parameter model matrix expression by parameter model matrix M direct inversion in example are as follows:
M=(α β ρ ε δ γ)T, (8)
According to gauss-newton method, first order Taylor is carried out near true value R and is unfolded to obtain:
J Δ M=Rk+1-Rk, (9)
Here k represents the number of iterations, Δ M=(Δ α Δ β Δ ρ Δ ε Δ δ Δ γ)TIt is the model updated, J is
The Jacobian matrix of HTI medium:
In (9) formula equal sign the right and left with multiplied by JT, it is available:
JTJ Δ M=JT(Rk+1-Rk), (11)
Enable H=JTJ brings the deformation of (11) formula into again and can obtain:
Δ M=H-1JT(Rk+1-Rk), (12)
Here for Hessian matrix, it is symmetrical and positive definite matrix to H.According to least square method, (9) formula is adjusted are as follows:
Δ M=(H+ λ I)-1JT(Rk+1-Rk), (13)
λ is a scalar, and I is unit matrix.It can establish an objective function according to above formula:
Q=| | (H+ λ I) Δ Μ-JT(Robs-RM)||2. (14)
Here RobsAnd RMRespectively observe reflection coefficient matrix and model reflection coefficient matrix.When objective function Q reaches
When satisfactory minimum value, iteration will stop.In this manner it is possible to obtain the model parameter closest to observation data.Gauss-ox
Pause algorithm, and the research of Sheen etc. illustrates it with very high resolution capability and convergence rate.
3, J matrix
It is different from parameter comparison inverting, the not clear expression formula of the partial derivative of Jacobian matrix J, J be it is big and sparse,
Precision will affect the result of inverting.Extension expression formula of the Jacobian matrix J in AVO are as follows:
Here reflection coefficient matrix and model parameter matrix is all n row matrix, and reflection coefficient is joined about elasticity in above formula
Several and anisotropic parameters single order local derviations, can be indicated with centered Finite Difference Methods:
With PP wave, PS1 wave, PS2 wave AVA joint inversion, then the J matrix in (10) formula is rewritten are as follows:
And H=JTJ rewrites (13) formula according to (17) formula are as follows:
E in above formulaPP、ePS1And ePS2It is PP wave respectively, the residual error of PS1 wave, PS2 wave AVO is all noises and its
The summation of his error.
4, constraint condition
Compared to traditional AVO inverting, combine the AVO inverting of PP, PS1, PS2 wave, using the data inversion of three kinds of waves
Six parameters out.Therefore, main to need to solve there are two problem: 1) to be compared with isotropism AVO inverting, AVO in anisotropy
Inverting is more likely to obtain local minimum solution, because the latter more three need the anisotropic parameters of inverting；2) as PP PS1
When the residual error of PS2 meets most reasonable value, the corresponding condition for stopping iteration.
The method that the present embodiment uses substep inverting, first to isotropic elasticity parameter: vertical, fast transverse wave speed and density
Inverting is carried out, second step continues inverting using the result of first step inverting as one of constraint condition.For first step inverting
As a result, up-and-down boundary is calculated as follows out in the minimum value a and maximum value b that obtain, its value is expanded 100% and is used as second step inverting
Constraint condition:
Inverting for anisotropic parameters is established initial each according to rock physics well logging information or velocity analysis
Its value is expanded 100% constraint condition as Anisotropic parameters inversion by anisotropy parameter model.
Due to truthful data and noise and other inevitable errors in joint inversion, increases three and stop
The only condition of iteration:
Ideally, the value of above formula is approximately 0.
5, trace gather is converted
Since the present embodiment is to carry out inverting to the AVA trace gather of Amplitudeversusangle rule, it is therefore desirable to prestack
Data is converted:
θ=arctan (Xp/ h), (21)
θ is incidence angle；XpFor the distance between shot point to reflection point；H is the depth of reflection point.It can similarly calculate
PS1 wave, PS2 wave correspond to incidence angle, then the amplitude value of corresponding data is assigned in different angle gathers.
6, time match
Since the ripple ratio longitudinal wave velocity of converted wave wants small, so hourage is longer, i.e., same reflection interface is come
Say, on the time section of three, PS1 wave, PS2 wave time ratio PP wave time it is big.Therefore, it is necessary to PP wave, PS1 wave and
PS2 wave profile carries out time match, carrys out adverse effect caused by release rate difference, the time of PS1 wave, PS2 wave is compressed to
The time of PP wave is upper:
tPS1(tPP)=tPP+u(tPP)
u(tPP) it is the conversion time difference, pass through each interval velocity ratio (VP-VS1)/VPIterative calculation, it is poor to obtain precise time, passes through
Practical to calculate, this method matching effect is good.
The process that the present embodiment provides a kind of fracture medium PP wave and division PS wave AVO joint inversion method is as follows:
S1, AVA trace gather extract
S1.1, PP wave AVA trace gather extract:
1.1.1 the depth of reflection point, and the corresponding incidence angle of different shot points) are found out according to the trace header information of seismic channel；
1.1.2) for needing the incidence angle of inverting, its amplitude value is extracted, new PP wave AVA trace gather is formed；
S1.2, PS1 wave AVA trace gather extract:
1.2.1 the depth of reflection point, and the corresponding incidence angle of different shot points) are found out according to the trace header information of seismic channel；
1.2.2) for needing the incidence angle of inverting, its amplitude value is extracted, new PS1 wave AVA trace gather is formed；
1.2.3 on when) joint velocity of longitudinal wave vertically travels PS1 wave data compression to PP wave, compressed PP wave is obtained
AVA trace gather, convenient for faults comparison between more waves:
tPS1(tPP)=tPP+u(tPP)
u(tPP) it is the conversion time difference, pass through speed ratio (VP-VS1)/VPIterative calculation, obtains the accurate time difference.
S1.3, PS2 wave AVA trace gather extract
S1.3.1, slow shear-wave S2 speed be not easy to obtain, be based onCalculate the speed of slow shear-wave；Further according to step
The information of PS1 wave in rapid S1.2 finds out the depth of reflection point, and the corresponding incidence angle of different shot points；
S1.3.2, the incidence angle for needing inverting, extract its amplitude value, form new PS2 wave AVA trace gather；
S1.3.3, similarly according to step S1.2.3 compressed PS2 wave AVA trace gather is obtained.
The foundation of S2, initial model
S2.1, high-resolution P- and S-wave velocity and density model are obtained by log data；
S2.2, according to rock physics well logging information or velocity analysis, establish initial anisotropic parameters model.
S3, inverting application
S3.1, anisotropic parameters model initial according to obtained in step S2 are calculated separately with formula (1) (2) (3)
PP wave, PS1 wave, PS2 wave reflection coefficient, anisotropic parameters are 0 in anisotropic parameters model initial at this time, not to each
Anisotropy parameter carries out inverting:
S3.2, using theoretical Ricker wavelet, the reflection system of the PP wave being calculated using step S3.1, PS1 wave, PS2 wave
Number makes the composite traces of PP wave, PS1 wave, PS2 wave respectively；
S3.3, combine PP wave, PS1 wave, PS2 wave calculating target function, and the essence of evaluation goal function using formula (14)
Degree；
S3.4, when objective function is unsatisfactory for required precision, then modify the model data of each layer, and return to step S3.1；
When objective function meets required precision, go to step S3.5；
S3.5, stop iteration；
S3.6, anisotropic parameters initial model is added in the inverting of step S3.1-S3.5；
S3.7, when objective function meets required precision, stop iteration, obtain at this time one reach well logging resolution ratio it is vertical,
Shear wave velocity and density and anisotropic parameters data volume.
The performance of the present embodiment inversion method is further illustrated below by way of test.
In order to test inversion technique provided by the invention, the present embodiment devises four layers of 1D HTI model,
Data are shown in Table 1.
Table 1
The reflection of PP wave, PS1 wave and PS2 wave is calculated based on the reflection coefficient approximate formula of R ü ger according to model data
Then coefficient obtains the AVA composite traces of PP wave, PS1 wave and PS2 wave with Ricker wavelet convolution, time compression is carried out to it
Match.Fig. 1 is the AVA trace gather of PP wave, PS1 wave and PS2 wave, and incidence angle is 4-24 °, is divided into 4 °, and (a) indicates PP wave, (b) indicates
PS1 wave (c) indicates PS2 wave.
In the test of this theoretical model, linear introductory die is made with the data of the first layer of parameter and the last layer
Type.The trace gather sparse for sampled point, inverting iteration speed are fast.In Fig. 2, dotted line indicates that true model, solid line indicate inverting
Model, dash line indicate initial model, and as can be seen from FIG. 2, inversion result and true model data are very close, anisotropy ginseng
Number ε(V)、δ(V)Inverting data essentially coincided with truthful data.
Inversion method provided in the present embodiment is suitable for crack HTI medium and isotropic medium, this is because instead
It penetrates in coefficient approximate formula, when anisotropic parameters are 0, as isotropic medium, even if initial model is linear, this reality
The inversion method for applying example can also obtain good effect.The multiwave AVO inversion method of HTI medium can be fracture medium stratum
Lithology and the effective informations such as fluid and fracture spacing are provided.
For those skilled in the art, it can be made various corresponding according to above technical solution and design
Change and modification, and all these change and modification should be construed as being included within the scope of protection of the claims of the present invention.
Claims (2)
1. a kind of fracture medium PP wave and division PS wave AVO joint inversion method, which comprises the steps of:
S1, PP wave AVA trace gather, PS1 wave AVA trace gather and PS2 wave AVA trace gather are extracted, and the time of PS1 wave, PS2 wave is compressed to
On the time of PP wave；
The foundation of S2, initial model: obtained by log data vertical speed, the vertical speed of fast transverse wave of initial longitudinal wave with
And density model, and according to rock physics well logging information or velocity analysis, establish initial anisotropic parameters model；
S3, inverting:
S3.1 carries out inverting to the vertical speed β and density p of isotropic elasticity parameter, that is, longitudinal wave vertical speed α, fast transverse wave:
S3.1.1, according to the model data of the vertical speed of initial longitudinal wave, the vertical speed of fast transverse wave and density model, use
The reflection coefficient approximate formula of R ü ger calculates separately the reflection coefficient of PP wave, PS1 wave, PS2 wave, at this time anisotropic parameters
ε(V)、δ(V), γ 0, not to anisotropic parameters carry out inverting；The reflection coefficient approximate formula of the R ü ger is as follows:
Wherein, Rpp、Rps1、Rps2The reflection coefficient of PP wave, PS1 wave, PS2 wave is respectively indicated, α indicates the vertical speed of longitudinal wave, β table
Show the vertical speed of fast transverse wave, ε(V)Refer to the difference of longitudinal wave horizontal velocity and vertical speed and the ratio of longitudinal wave vertical speed, indicates
The anisotropic degree of longitudinal wave；δ(V)Indicate longitudinal wave and shear wave anisotropy relative size；The anisotropy of γ expression shear wave velocity
Development degree；Z indicates that vertical p-wave impedance Z=ρ α, G indicate shear modulus G=ρ β of the wave of corresponding vertical transmission2, i is vertical
Wave incidence angle, j are fast transverse wave angle of reflection, and θ is incident orientation angle, and subscript 1 represents the physical parameter of top dielectric, and subscript 2 represents
The parameter of layer dielectric；
S3.1.2, using theoretical Ricker wavelet, calculated PP wave, PS1 wave, PS2 wave reflection coefficient make respectively
PP wave, PS1 wave, PS2 wave composite traces；
S3.1.3, joint PP wave, PS1 wave, PS2 wave calculating target function are as follows:
Q=| | (H+ λ I) Δ Μ-JT(Robs-RM)||2；
Wherein, M=(α β ρ ε δ γ)TSingle order is carried out near true value R according to gauss-newton method for parameter model matrix
Taylor expansion obtains J Δ M=Rk+1-Rk；K represents the number of iterations, Δ M=(Δ α Δ β Δ ρ Δ ε Δ δ Δ γ)TIt is to update
Model, J be HTI medium Jacobian matrix；RobsAnd RMRespectively observe reflection coefficient matrix and model reflection coefficient square
Battle array；λ is a scalar, and I is unit matrix；
S3.1.4, when objective function is unsatisfactory for required precision, then modify the vertical speed of longitudinal wave, the vertical speed of fast transverse wave and
The model data of density model, and the reflecting system of return step S3.1.1 calculates step again, carries out inverting again；Work as target
Function meets required precision, stops iteration；Finally obtain the vertical speed α, fast for reaching the inversion result longitudinal wave of well logging resolution ratio
The vertical speed β and density p of shear wave；
S3.2, using the inversion result of step S3.1 as constraint condition, to anisotropic parameters carry out inverting:
S3.2.1, using the model data of initial anisotropic parameters model, distinguished with the reflection coefficient approximate formula of R ü ger
The reflection coefficient for calculating PP wave, PS1 wave, PS2 wave, in the reflection coefficient approximate formula of the R ü ger:
S3.2.2, using theoretical Ricker wavelet, calculated PP wave, PS1 wave, PS2 wave reflection coefficient make respectively
PP wave, PS1 wave, PS2 wave composite traces；
S3.2.3, combine PP wave, PS1 wave, PS2 wave calculating target function according to the method for step S3.1.3；
S3.2.4, when objective function is unsatisfactory for required precision, then modify the vertical speed of longitudinal wave, the vertical speed of fast transverse wave and
The model data of density model, and the reflecting system of return step S3.2.1 calculates step again, carries out inverting again；Work as target
Function meets required precision, stops iteration；Finally obtain the inversion result anisotropic parameters ε for reaching well logging resolution ratio(V)、
δ(V),γ；
Wherein,
In the objective function, extension expression formula of the Jacobian matrix J in AVO are as follows:
Single order local derviation of the reflection coefficient about elastic parameter and anisotropic parameters, is indicated with centered Finite Difference Methods:
With PP wave, PS1 wave, PS2 wave AVA joint inversion, then Jacobian matrix J are as follows:
And H=JTJ then has:
Wherein, ePP、ePS1And ePS2It is PP wave, the residual error of PS1 wave, PS2 wave AVO, θ respectively1…θiRepresent each of longitudinal wave
A incidence angle, R1…RiRepresent the reflection coefficient of PP wave, PS1 wave corresponding to each angle, PS2 wave.
2. fracture medium PP wave according to claim 1 and division PS wave AVO joint inversion method, which is characterized in that
The detailed process of step S1 are as follows:
S1.1, PP wave AVA trace gather extract；
1.1.1 the depth of reflection point, and the corresponding incidence angle of different shot points) are found out according to the trace header information of seismic channel；
1.1.2) for needing the incidence angle of inverting, its amplitude value is extracted, new PP wave AVA trace gather is formed；
S1.2, PS1 wave AVA trace gather extract:
1.2.1 the depth of reflection point, and the corresponding incidence angle of different shot points) are found out according to the trace header information of seismic channel；
1.2.2) for needing the incidence angle of inverting, its amplitude value is extracted, new PS1 wave AVA trace gather is formed；
1.2.3) time of PS1 wave is compressed on the vertical hourage of PP wave by joint velocity of longitudinal wave, obtains compressed PS1
Wave AVA trace gather:
tPS1(tPP)=tPP+u(tPP)；
tppIndicate the vertical hourage of PP wave, tps1(tpp) indicate PS1 wave time；U (tPP) is the conversion time difference；
S1.3, PS2 wave AVA trace gather extract
S1.3.1, it is based onCalculate the speed of slow shear-wave；Wherein, γ indicates that the anisotropy of shear wave velocity develops journey
Degree, Vs2Indicate the speed of slow shear-wave, Vs1Indicate the speed of fast transverse wave；It is found out instead further according to the information of the PS1 wave in step S1.2
The depth of exit point, and the corresponding incidence angle of different shot points；
S1.3.2, the incidence angle for needing inverting, extract its amplitude value, form new PS2 wave AVA trace gather；
The time of PS2 wave is compressed on the vertical hourage of PP wave by S1.3.3, joint velocity of longitudinal wave, obtains compressed
PS2 wave AVA trace gather:
tps2(tpp)=tpp+u(tpp)；
tppIndicate the vertical hourage of PP wave, tps2(tpp) indicate PS2 wave time；U (tPP) is the conversion time difference.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810894632.5A CN108828660B (en) | 2018-08-08 | 2018-08-08 | A kind of fracture medium PP wave and division PS wave AVO joint inversion method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810894632.5A CN108828660B (en) | 2018-08-08 | 2018-08-08 | A kind of fracture medium PP wave and division PS wave AVO joint inversion method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108828660A CN108828660A (en) | 2018-11-16 |
| CN108828660B true CN108828660B (en) | 2019-10-22 |
Family
ID=64152697
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810894632.5A Active CN108828660B (en) | 2018-08-08 | 2018-08-08 | A kind of fracture medium PP wave and division PS wave AVO joint inversion method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108828660B (en) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110187389B (en) * | 2019-06-24 | 2020-07-24 | 中国地质大学（北京） | AVA inversion method based on thin layer reflection theory |
| CN112198549B (en) * | 2019-07-08 | 2024-05-28 | 中国石油天然气集团有限公司 | Pre-stack crack determination method and system based on seismic forward modeling board |
| CN111551990B (en) * | 2020-04-07 | 2023-05-23 | 西安科技大学 | Acquisition system for HTI type coal seam seismic wave reflection coefficient |
| CN112230276B (en) * | 2020-09-30 | 2023-07-07 | 兰州石化职业技术学院 | Crack type tight reservoir fluid identification method, system, identification instrument, medium and application |
Family Cites Families (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6785612B1 (en) * | 2003-05-29 | 2004-08-31 | Pgs Americas, Inc. | Seismic velocity update for anisotropic depth migration |
| CN102053277A (en) * | 2009-10-30 | 2011-05-11 | 中国石油化工股份有限公司 | Method for detecting reservoir fissure development direction by utilizing seismic data |
| CN103163554A (en) * | 2013-02-04 | 2013-06-19 | 西安交通大学 | Self-adapting wave form retrieval method through utilization of zero offset vertical seismic profile (VSP) data to estimate speed and Q value |
| WO2016027156A1 (en) * | 2014-08-19 | 2016-02-25 | Cgg Services Sa | Joint inversion of compressional and shear seismic data in native time domains |
| CN104965224B (en) * | 2015-06-03 | 2017-10-27 | 北京多分量地震技术研究院 | PP ripples, which are carried out, with average incident angle gathers combines AVO inversion methods with PS ripples |
| CN106033127B (en) * | 2016-06-29 | 2018-04-06 | 中国石油化工股份有限公司 | Crustal stress azimuthal seismic Forecasting Methodology based on shear wave velocity rate of change |
-
2018
- 2018-08-08 CN CN201810894632.5A patent/CN108828660B/en active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN108828660A (en) | 2018-11-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108828660B (en) | A kind of fracture medium PP wave and division PS wave AVO joint inversion method | |
| EP3076205B1 (en) | Method for survey data processing compensating for visco-acoustic effects in tilted transverse isotropy reverse time migration | |
| Binder et al. | Modeling the seismic response of individual hydraulic fracturing stages observed in a time-lapse distributed acoustic sensing vertical seismic profiling survey | |
| EP3028071B1 (en) | Method and device for the generation and application of anisotropic elastic parameters in horizontal transverse isotropic (hti) media | |
| US10520619B2 (en) | FWI model domain angle stacks with amplitude preservation | |
| CN104122585A (en) | Seismic forward modeling method based on elastic wave field vector decomposition and low-rank decomposition | |
| CN109669212B (en) | Seismic data processing method, stratum quality factor estimation method and device | |
| CN108345031A (en) | A kind of elastic fluid active source and passive source, which are mixed, adopts seismic data full waveform inversion method | |
| CA2386568A1 (en) | System for estimating azimuthal variations in seismic data | |
| WO2012139082A1 (en) | Event selection in the image domain | |
| CN110780351B (en) | Longitudinal wave and converted wave prestack joint inversion method and system | |
| CN111722284B (en) | Method for establishing speed depth model based on gather data | |
| CN104237945A (en) | Seismic data self-adaptive high-resolution processing method | |
| CN113031068B (en) | Reflection coefficient accurate base tracking prestack seismic inversion method | |
| CN105093278A (en) | Extraction method for full waveform inversion gradient operator based on excitation main energy optimization algorism | |
| CN109188511A (en) | A kind of thin sand-mud interbed medium multi-wave AVO joint inversion method | |
| US10955576B2 (en) | Full waveform inversion of vertical seismic profile data for anisotropic velocities using pseudo-acoustic wave equations | |
| CN106501872B (en) | A kind of calculation method and device of fracture reservoir ground stress characteristics | |
| CN111208564A (en) | Depth domain horizon calibration method and device | |
| Mizuno et al. | Anisotropic velocity model inversion for imaging the microseismic cloud | |
| US20190346581A1 (en) | Methods for determining transversely isotropic-elastic constants from borehole sonic velocities in strongly transversely-isotropic formations | |
| CN111025388B (en) | Multi-wave combined prestack waveform inversion method | |
| Kazei et al. | Amplitude-based DAS logging: Turning DAS VSP amplitudes into subsurface elastic properties | |
| Li et al. | Active-source Rayleigh wave dispersion by the Aki spectral formulation | |
| WO2015124960A1 (en) | Systems and methods for improved inversion analysis of 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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |