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
- waves
- longitudinal
- ava
- inversion
- 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 40
- 239000011159 matrix material Substances 0.000 claims description 29
- 238000004364 calculation method Methods 0.000 claims description 9
- 230000014509 gene expression Effects 0.000 claims description 7
- 238000011161 development Methods 0.000 claims description 5
- 239000011435 rock Substances 0.000 claims description 5
- 238000000605 extraction Methods 0.000 claims description 4
- 230000000704 physical effect Effects 0.000 claims description 3
- 230000001902 propagating effect Effects 0.000 claims description 3
- 230000007704 transition Effects 0.000 claims 2
- 230000006835 compression Effects 0.000 abstract description 5
- 238000007906 compression Methods 0.000 abstract description 5
- 230000000694 effects Effects 0.000 abstract description 5
- 239000012530 fluid Substances 0.000 abstract description 2
- 238000012360 testing method Methods 0.000 description 5
- 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
- 238000005070 sampling 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
- 230000002411 adverse Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 230000009466 transformation 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 invention relates to the technical field of seismic exploration, in particular to a joint inversion method of a PP wave and a split PS wave AVO (amplitude versus offset) of a crack medium.
Background
The AVO (amplitude variation with offset)/AVA (amplitude variation with angle) inversion technique is based on elastic wave theory, and utilizes the change rule of reflection amplitude of prestack seismic gathers along with offset (or incident angle) to estimate and predict physical parameters and rock characteristics of reservoir.
The crack-induced HTI medium is a theoretical equivalent medium to inlay a set of vertically aligned cracks in an isotropic medium, which is an anisotropic medium with a transverse horizontal axis of symmetry. In HTI media, AVO inversion requires more parameters because the anisotropy of the media causes the reflection coefficient to change, and the AVO inversion also needs to take into account the anisotropy parameters of the media. For the AVO inversion research of the anisotropic medium, due to the large number of parameters and complexity, the multi-solution property of the conventional AVO inversion method cannot be solved, and further, a good prediction effect cannot be obtained. Lee et al propose a method that more easily solves the problem of anisotropic media with local minimum solutions; auger et al, Verie and Landro propose joint inversion of PP waves and PS waves to solve the problem of multi-solution, and obtain better results in isotropic media; fatti et al accurately invert longitudinal and transverse wave impedances by using a weighted stacking method; nadri and Hartly invert elastic parameters and anisotropic parameters by using a conjugate iterative algorithm; lu et al jointly invert 5 items of anisotropy parameters and elasticity parameters for VTI medium multi-wave.
Disclosure of Invention
Aiming at the defects of the prior art, the invention aims to provide an AVO (amplitude versus offset) joint inversion method for a PP wave and a split PS wave of a fracture medium, which is characterized in that time compression matching is carried out on an AVA angle gather before folding, and then based on an R ü ger reflection coefficient approximate formula about an HTI medium, a Gauss Newton iterative algorithm is applied to carry out joint inversion on the PP wave, the PS1 wave and the PS2 wave, so that the results of anisotropic parameters and elastic parameters of the fracture medium are obtained.
In order to achieve the purpose, the invention adopts the following technical scheme:
a fracture medium PP wave and split PS wave AVO joint inversion method comprises the following steps:
s1, extracting a PP wave AVA track set, a PS1 wave AVA track set and a PS2 wave AVA track set, and compressing the time of PS1 waves and PS2 waves to the time of the PP waves;
s2, establishing an initial model: obtaining a vertical speed of an initial longitudinal wave, a vertical speed of a fast transverse wave and a density model through logging data, and establishing an initial anisotropic parameter model according to rock physical logging information or speed analysis;
s3, inversion:
s3.1, inverting the isotropic elastic parameters, namely the vertical speed alpha of longitudinal waves, the vertical speed beta of fast transverse waves and the density rho:
s3.1.1, calculating the reflection coefficients of PP wave, PS1 wave and PS2 wave respectively by R ü ger reflection coefficient approximate formula according to the initial vertical velocity of longitudinal wave, the vertical velocity of fast transverse wave and the model data of density model, wherein the anisotropic parameter epsilon is(V)、δ(V)Gamma is 0, no inversion is performed on the anisotropy parameters, the inverse of R ü gerThe equation for the approximation of the coefficient of radiation is as follows:
wherein R ispp、Rps1、Rps2Respectively represents the reflection coefficients of PP wave, PS1 wave and PS2 wave, alpha represents the vertical velocity of longitudinal wave, beta represents the vertical velocity of fast transverse wave, epsilon(V)The ratio of the difference value of the horizontal velocity and the vertical velocity of the longitudinal wave to the vertical velocity of the longitudinal wave represents the anisotropy degree of the longitudinal wave; delta(V)Representing the relative magnitude of the anisotropy of longitudinal waves and transverse waves; γ represents the degree of anisotropic development of the shear wave velocity; z denotes the vertical longitudinal wave impedance Z ═ ρ α, and G denotes the shear modulus G ═ ρ β of the corresponding vertically propagating wave2I is a longitudinal wave incident angle, j is a fast transverse wave reflection angle, theta is an incident azimuth angle, subscript 1 represents a physical property parameter of an upper medium, and subscript 2 represents a parameter of a lower medium;
s3.1.2, respectively making synthetic records of PP waves, PS1 waves and PS2 waves by using the reflection coefficients of the PP waves, the PS1 waves and the PS2 waves obtained by calculation by adopting theoretical Rake wavelets;
s3.1.3, calculating an objective function by combining the PP wave, the PS1 wave and the PS2 wave as follows:
Q＝||(H+λI)ΔΜ-JT(Robs-RM)||2；
where M ═ α β ρ ∈ δ γ)TPerforming first-order Taylor expansion near a true value R according to a Gauss-Newton method to obtain a parameter model matrixk+1-Rk(ii) a k represents the number of iterations, Δ M ═ Δ α Δ β Δ ρ Δ ∈ Δ γ)TIs an updated model, J is the Jacobian matrix of the HTI medium; robsAnd RMRespectively is the observed reflection coefficientA matrix and a model reflection coefficient matrix;
s3.1.4, when the target function does not meet the precision requirement, modifying the vertical velocity of longitudinal waves, the vertical velocity of fast transverse waves and the model data of the density model, and returning to the reflecting system calculation step of step S3.1.1 again to perform inversion again; when the target function meets the precision requirement, stopping iteration; finally, the vertical speed alpha of longitudinal waves, the vertical speed beta of fast transverse waves and the density rho of inversion results reaching the logging resolution are obtained;
s3.2, taking the inversion result of the step S3.1 as a constraint condition, and inverting the anisotropic parameters:
s3.2.1, respectively calculating the reflection coefficients of the PP wave, the PS1 wave and the PS2 wave by using the model data of the initial anisotropic parameter model and a reflection coefficient approximation formula of R ü ger;
s3.2.2, respectively making synthetic records of PP waves, PS1 waves and PS2 waves by using the reflection coefficients of the PP waves, the PS1 waves and the PS2 waves obtained by calculation by adopting theoretical Rake wavelets;
s3.2.3, calculating an objective function by combining the PP wave, the PS1 wave and the PS2 wave according to the method of the step S3.1.3;
s3.2.4, when the target function does not meet the precision requirement, modifying the vertical velocity of longitudinal waves, the vertical velocity of fast transverse waves and the model data of the density model, and returning to the reflecting system calculation step of step S3.2.1 again to perform inversion again; when the target function meets the precision requirement, stopping iteration; finally obtaining the anisotropy parameter epsilon of the inversion result reaching the logging resolution(V)、δ(V)、γ。
Further, the specific process of step S1 is:
s1.1, extracting a PP wave AVA gather;
1.1.1) calculating the depth of a reflection point and the corresponding incident angle of different shot points according to the track head information of the seismic channel;
1.1.2) extracting amplitude values of incident angles needing inversion to form a new PP wave AVA gather;
s1.2, extracting an AVA trace set of PS1 waves:
1.2.1) calculating the depth of a reflection point and the corresponding incident angle of different shot points according to the track head information of the seismic channel;
1.2.2) extracting amplitude values of incidence angles needing inversion to form a new PS1 wave AVA gather;
1.2.3) compressing the time of the PS1 wave to the vertical travel time of the PP wave by combining the longitudinal wave velocity to obtain a compressed PS1 wave AVA gather:
tPS1(tPP)＝tPP+u(tPP)；
tpprepresenting the vertical travel time, t, of the PP-waveps1(tpp) Time representing PS1 waves; u (t)PP) Converting the time difference;
s1.3, PS2 wave AVA gather extraction
S1.3.1, based onCalculating the speed of the slow transverse wave; wherein γ represents the degree of anisotropic development of the transverse wave velocity, Vs2Indicating the speed, V, of the slow transverse waves1Representing the velocity of the fast transverse wave; then, the depth of a reflection point and the incident angle corresponding to different shot points are obtained according to the information of the PS1 wave in the step S1.2;
s1.3.2, extracting amplitude values of the incident angles to be inverted to form a new PS2 wave AVA gather;
s1.3.3, compressing the time of the PS2 wave to the vertical travel time of the PP wave by combining with the longitudinal wave velocity to obtain a compressed PS2 wave AVA gather:
tps2(tpp)＝tpp+u(tpp)；
tpprepresenting the vertical travel time, t, of the PP-waveps2(tpp) Time representing PS2 waves; u (t)PP) To convert the moveout.
Further, the reflection coefficient of R ü ger is approximated by the formula:
further, in the objective function, an extended expression of the jacobian matrix J in the AVO is:
the first-order partial derivatives of the reflection coefficients with respect to the elastic and anisotropic parameters are expressed by the central difference method:
by applying the joint inversion of PP wave, PS1 wave and PS2 wave AVA, the Jacobian matrix J is as follows:
and H ═ JTJ, then there are:
wherein e isPP、ePS1And ePS2The residual errors of the PP wave, PS1 wave, and PS2 wave AVO, respectively, are the sum of all noise and other errors.
The invention has the beneficial effects that:
according to the method, time compression matching is carried out on the pre-stack AVA angle gather, then based on an R ü ger reflection coefficient approximate formula about an HTI medium, Gaussian Newton iteration algorithm is applied, PP wave, PS1 wave and PS2 wave combined inversion is carried out, so that anisotropic parameter and elastic parameter results of the fracture medium are obtained, and the method is suitable for the fracture HTI medium and isotropic medium.
Drawings
FIG. 1 is a schematic diagram of AVA gathers of a PP wave, a PS1 wave, and a PS2 wave under test in accordance with an embodiment of the present invention;
FIG. 2 is a schematic diagram of inversion results in a test according to an embodiment of the present invention.
Detailed Description
The present invention will be further described with reference to the accompanying drawings, and it should be noted that the present embodiment provides detailed embodiments and specific operation procedures based on the premise of the technical solution, but the protection scope of the present invention is not limited to the present embodiment.
The principles of the method provided by the present embodiment will be described and illustrated in detail below.
1. Coefficient of reflection
In the crack HTI medium, transverse waves are split into fast transverse waves with the vibration direction parallel to the crack surface and slow transverse waves with the vibration direction perpendicular to the crack surface, so that the AVO inversion is facilitated by using the display expressions of the reflection coefficients of PP, PS1 and PS2 waves for the crack HTI medium, in the embodiment, the HTI medium is inverted by using the approximate formula of R ü ger:
where α represents the vertical velocity of longitudinal waves, β represents the vertical velocity of fast transverse waves, and ε(V)The ratio of the difference value of the horizontal velocity and the vertical velocity of the longitudinal wave to the vertical velocity of the longitudinal wave represents the anisotropy degree of the longitudinal wave; delta(V)Representing the relative magnitude of the anisotropy of longitudinal waves and transverse waves; γ represents the degree of anisotropic development of the shear wave velocity. Where Z denotes the vertical longitudinal wave impedance Z ═ ρ α, and G denotes the shear modulus G ═ ρ β of the corresponding vertically propagating wave2I is the incident angle of longitudinal wave, j is the reflection angle of fast transverse wave, and theta isThe subscripts PP, PS1, and PS2 indicate plane wave types in which the incident wave is a longitudinal wave and the reflected wave is a longitudinal wave, a fast transverse wave, and a slow transverse wave, respectively. Subscript 1 represents the physical property parameters of the upper medium, subscript 2 represents the parameters of the lower medium, and in the above formula:
according to the reflection coefficient, the time window of n sampling points is used for inversion, and the matrix expression of the reflection coefficient is as follows:
theta here1…θiRepresents each incident angle, R, of the longitudinal wave1…RiThe reflection coefficients of the PP wave, the PS1 wave and the PS2 wave corresponding to each angle are represented. Thus, PP wave, PS1 wave, PS2 wave AVA seismic data are:
S＝WR， (6)
w is the wavelet matrix, which is described as:
for the AVA inversion of the HTI medium, a Gaussian-Newton method is used for solving the nonlinear least square problem. Given observation B ═ B1…bm) And n sets of model parameters X ═ X (X)1…xn) Initially guessing model data X(0)As a first iteration parameter, the gauss-newton algorithm can make the model parameter approach the true value by iteratively calculating the minimum variance Q.
2. Model iteration
The model parameters of elastic contrast are used in the R ü ger (2002) approximation formulaHowever, in this embodiment, the attributes of each layer are directly inverted through the parametric model matrix M, so the parametric model matrix MThe array expression is:
M＝(α β ρ ε δ γ)T, (8)
according to the Gauss-Newton method, performing first-order Taylor expansion near the true value R to obtain:
JΔM＝Rk+1-Rk， (9)
where k denotes the number of iterations, Δ M ═ Δ α Δ β Δ ρ Δ ∈ Δ γ)TIs an updated model, J is the jacobian matrix of the HTI medium:
multiplying the left side and the right side of the equation (9) by JTIt is possible to obtain:
JTJΔM＝JT(Rk+1-Rk), (11)
let H equal JTThe transformation of J into formula (11) can be obtained:
ΔM＝H-1JT(Rk+1-Rk), (12)
here H is the hessian matrix, which is a symmetric and positive definite matrix. According to the least square method, the formula (9) is adjusted as follows:
ΔM＝(H+λI)-1JT(Rk+1-Rk), (13)
λ is a scalar and I is an identity matrix. An objective function can be established according to the above equation:
Q＝||(H+λI)ΔΜ-JT(Robs-RM)||2. (14)
r in this caseobsAnd RMRespectively an observed reflection coefficient matrix and a model reflection coefficient matrix. When the objective function Q reaches a satisfactory minimum, the iteration will stop. In this way, the model parameters closest to the observed data can be obtained. The research of Gauss-Newton algorithm, Sheen and the like shows that the method has high resolution capability and convergence rate.
3. J matrix
Unlike parametric contrast inversion, the partial derivatives of the jacobian matrix J have no explicit expression, J is large and sparse, and its accuracy will affect the inversion result. The extended expression of the Jacobian matrix J in AVO is as follows:
the reflection coefficient matrix and the model parameter matrix are n-row matrices, and the first-order partial derivatives of the reflection coefficients with respect to the elastic parameters and the anisotropic parameters in the above formula can be expressed by a central difference method:
by applying the joint inversion of PP wave, PS1 wave and PS2 wave AVA, the J matrix in the formula (10) is rewritten as:
and H ═ JTJ, according to equation (17), rewrite equation (13) as:
e in the above formulaPP、ePS1And ePS2The residual errors of the PP wave, PS1 wave, and PS2 wave AVO, respectively, are the sum of all noise and other errors.
4. Constraint conditions
Compared with the traditional AVO inversion, the method combines the AVO inversion of PP, PS1 and PS2 waves, and adopts the data inversion of three waves to obtain six parameters. Thus, there are two main problems to be solved: 1) compared with isotropic AVO inversion, in anisotropy, AVO inversion is more prone to obtain local minimum solutions, because the latter has three more anisotropy parameters to be inverted; 2) when the residual error of PP PS1PS2 satisfies the most reasonable value, the corresponding condition for stopping iteration is satisfied.
In this embodiment, a step-by-step inversion method is adopted, and first, for isotropic elastic parameters: and inverting the longitudinal and fast transverse wave speeds and the density, and continuing to invert by taking the result of the first-step inversion as one of the constraint conditions in the second step. For the result of the first-step inversion, the obtained minimum value a and the maximum value b calculate the upper and lower boundaries according to the following formula, and the values are expanded by 100% to be used as the constraint conditions of the second-step inversion:
for the inversion of the anisotropic parameters, an initial anisotropic parameter model is established according to the rock physical logging information or velocity analysis, and the value of the model is expanded by 100% to be used as a constraint condition for the inversion of the anisotropic parameters.
Three conditions for stopping the iteration are added in the joint inversion due to real data and noise and other unavoidable errors:
ideally, the value of the above equation is approximately 0.
5. Gather conversion
Since the present embodiment inverts the AVA gather whose amplitude varies with the incident angle, it is necessary to convert the prestack data:
θ＝arctan(Xp/h)， (21)
theta is an incident angle; xpThe distance between the shot point and the reflection point is taken as the distance; h is the depth of the reflection point. Similarly, the corresponding incident angles of the PS1 wave and the PS2 wave can be calculated, and the amplitude values of the corresponding data are distributed to different angle channel sets.
6. Time matching
The travel time is longer because the wave speed of the converted wave is lower than that of the longitudinal wave, namely, the time of the PS1 wave and the PS2 wave is longer than that of the PP wave on the time section of the three waves in the same reflection interface. Therefore, the PP wave, the PS1 wave and the PS2 wave profiles must be time-matched to eliminate the adverse effects caused by the speed difference, and the time of the PS1 wave and the PS2 wave is compressed to the time of the PP wave:
tPS1(tPP)＝tPP+u(tPP)
u(tPP) For converting time differences, the velocity ratio (V) of the layers is passedP-VS1)/VPIterative computation is carried out to obtain accurate time difference, and the matching effect of the method is good through actual computation.
The flow of the crack medium PP wave and split PS wave AVO joint inversion method provided by the embodiment is as follows:
s1, AVA gather extraction
S1.1, extracting a PP wave AVA gather:
1.1.1) calculating the depth of a reflection point and the corresponding incident angle of different shot points according to the track head information of the seismic channel;
1.1.2) extracting amplitude values of incident angles needing inversion to form a new PP wave AVA gather;
s1.2, extracting an AVA trace set of PS1 waves:
1.2.1) calculating the depth of a reflection point and the corresponding incident angle of different shot points according to the track head information of the seismic channel;
1.2.2) extracting amplitude values of incidence angles needing inversion to form a new PS1 wave AVA gather;
1.2.3) combine the longitudinal wave speed to compress the PS1 wave data to the PP wave vertical travel time, obtain the PP wave AVA trace set after compression, the level contrast between the multiwave of being convenient for:
tPS1(tPP)＝tPP+u(tPP)
u(tPP) For converting time differences, the velocity ratio (V) is passedP-VS1)/VPAnd (5) performing iterative computation to obtain the accurate time difference.
S1.3, PS2 wave AVA gather extraction
S1.3.1, the speed of the slow transverse wave S2 is not easy to obtain, based onCalculating the speed of the slow transverse wave; then, the depth of a reflection point and the incident angle corresponding to different shot points are obtained according to the information of the PS1 wave in the step S1.2;
s1.3.2, extracting amplitude values of the incident angles to be inverted to form a new PS2 wave AVA gather;
s1.3.3, obtaining the compressed PS2 wave AVA gather according to the step S1.2.3.
S2, establishing initial model
S2.1, obtaining high-resolution longitudinal and transverse wave speeds and density models through logging data;
and S2.2, establishing an initial anisotropic parameter model according to the rock physical logging information or speed analysis.
S3 application of inversion
S3.1, respectively calculating the reflection coefficients of the PP wave, the PS1 wave and the PS2 wave by using the formulas (1), (2) and (3) according to the initial anisotropic parameter model obtained in the step S2, wherein the anisotropic parameter in the initial anisotropic parameter model is 0, and the anisotropy parameter is not inverted:
s3.2, respectively making synthetic records of the PP wave, the PS1 wave and the PS2 wave by using the reflection coefficients of the PP wave, the PS1 wave and the PS2 wave which are obtained by calculation in the step S3.1 by adopting theoretical Rake wavelets;
s3.3, calculating a target function by combining the PP wave, the PS1 wave and the PS2 wave according to a formula (14), and evaluating the precision of the target function;
s3.4, when the target function does not meet the precision requirement, modifying the model data of each layer, and returning to the step S3.1; when the target function meets the precision requirement, skipping to the step S3.5;
s3.5, stopping iteration;
s3.6, adding the initial model of the anisotropic parameters into the inversion of the steps S3.1-S3.5;
and S3.7, stopping iteration when the target function meets the precision requirement, and obtaining a longitudinal wave speed, a transverse wave speed, a density and an anisotropic parameter data volume which reach the logging resolution at the moment.
The performance of the inversion method of the present embodiment is further illustrated by the following tests.
In order to test the inversion technique provided by the present invention, a four-layer 1D HTI model is designed in this embodiment, and the data is shown in table 1.
TABLE 1
According to model data, reflection coefficients of a PP wave, a PS1 wave and a PS2 wave are calculated based on a reflection coefficient approximation formula of R ü ger, and then convolution with a Rake wavelet is carried out to obtain AVA synthetic records of the PP wave, the PS1 wave and the PS2 wave, and time compression matching is carried out on the AVA synthetic records.
In the theoretical model test, a linear initial model is made by using the data of the first layer and the last layer of the parameters. And for a gather with sparse sampling points, the inversion iteration speed is high. In fig. 2, the dashed line represents the real model, the solid line represents the inverse model, and the dashed line represents the initial model, and it can be known from fig. 2 that the inverse result is very close to the real model data, and the anisotropy parameter epsilon(V)、δ(V)Substantially coinciding with the real data.
The inversion method provided by the embodiment is suitable for a crack HTI medium and an isotropic medium, and is because in the reflection coefficient approximation formula, when the anisotropy parameter is 0, the reflection coefficient is the isotropic medium, and even if the initial model is linear, the inversion method of the embodiment can also obtain a good effect. The multi-wave AVO inversion method of the HTI medium can provide effective information such as lithology, fluid, fracture density and the like for the fracture medium stratum.
Various changes and modifications can be made by those skilled in the art based on the above technical solutions and concepts, and all such changes and modifications should be included in the scope of the present invention.
Claims (2)
1. A joint inversion method of a crack medium PP wave and a split PS wave AVO is characterized by comprising the following steps:
s1, extracting a PP wave AVA track set, a PS1 wave AVA track set and a PS2 wave AVA track set, and compressing the time of PS1 waves and PS2 waves to the time of the PP waves;
s2, establishing an initial model: obtaining a vertical speed of an initial longitudinal wave, a vertical speed of a fast transverse wave and a density model through logging data, and establishing an initial anisotropic parameter model according to rock physical logging information or speed analysis;
s3, inversion:
s3.1, inverting the isotropic elastic parameters, namely the vertical speed alpha of longitudinal waves, the vertical speed beta of fast transverse waves and the density rho:
s3.1.1, calculating the reflection coefficients of PP wave, PS1 wave and PS2 wave respectively by R ü ger reflection coefficient approximate formula according to the initial vertical velocity of longitudinal wave, the vertical velocity of fast transverse wave and the model data of density model, wherein the anisotropic parameter epsilon is(V)、δ(V)Gamma is 0, and the anisotropy parameter is not inverted, and the reflection coefficient approximation formula of R ü ger is as follows:
wherein R ispp、Rps1、Rps2Respectively represents the reflection coefficients of PP wave, PS1 wave and PS2 wave, alpha represents the vertical velocity of longitudinal wave, beta represents the vertical velocity of fast transverse wave, epsilon(V)The ratio of the difference value of the horizontal velocity and the vertical velocity of the longitudinal wave to the vertical velocity of the longitudinal wave represents the anisotropy degree of the longitudinal wave; delta(V)Representing the relative magnitude of the anisotropy of longitudinal waves and transverse waves; γ represents the degree of anisotropic development of the shear wave velocity; z denotes the vertical longitudinal wave impedance Z ═ ρ α, and G denotes the shear modulus G ═ ρ β of the corresponding vertically propagating wave2I is the incident angle of longitudinal wave, j is the reflection angle of fast transverse wave, and theta is the incident angleThe azimuth angle, subscript 1 represents the physical property parameter of the upper medium, and subscript 2 represents the parameter of the lower medium;
s3.1.2, respectively making synthetic records of PP waves, PS1 waves and PS2 waves by using the reflection coefficients of the PP waves, the PS1 waves and the PS2 waves obtained by calculation by adopting theoretical Rake wavelets;
s3.1.3, calculating an objective function by combining the PP wave, the PS1 wave and the PS2 wave as follows:
Q＝||(H+λI)ΔΜ-JT(Robs-RM)||2；
where M ═ α β ρ ∈ δ γ)TPerforming first-order Taylor expansion near a true value R according to a Gauss-Newton method to obtain a parameter model matrixk+1-Rk(ii) a k represents the number of iterations, Δ M ═ Δ α Δ β Δ ρ Δ ∈ Δ γ)TIs an updated model, J is the Jacobian matrix of the HTI medium; robsAnd RMRespectively an observation reflection coefficient matrix and a model reflection coefficient matrix; λ is a scalar, and I is an identity matrix;
s3.1.4, when the target function does not meet the precision requirement, modifying the vertical velocity of longitudinal waves, the vertical velocity of fast transverse waves and the model data of the density model, and returning to the reflecting system calculation step of step S3.1.1 again to perform inversion again; when the target function meets the precision requirement, stopping iteration; finally, the vertical speed alpha of longitudinal waves, the vertical speed beta of fast transverse waves and the density rho of inversion results reaching the logging resolution are obtained;
s3.2, taking the inversion result of the step S3.1 as a constraint condition, and inverting the anisotropic parameters:
s3.2.1, respectively calculating the reflection coefficients of the PP wave, the PS1 wave and the PS2 wave by using R ü ger reflection coefficient approximate formulas by using model data of an initial anisotropic parameter model, wherein in the R ü ger reflection coefficient approximate formulas:
s3.2.2, respectively making synthetic records of PP waves, PS1 waves and PS2 waves by using the reflection coefficients of the PP waves, the PS1 waves and the PS2 waves obtained by calculation by adopting theoretical Rake wavelets;
s3.2.3, calculating an objective function by combining the PP wave, the PS1 wave and the PS2 wave according to the method of the step S3.1.3;
s3.2.4, when the target function does not meet the precision requirement, modifying the vertical velocity of longitudinal waves, the vertical velocity of fast transverse waves and the model data of the density model, and returning to the reflecting system calculation step of step S3.2.1 again to perform inversion again; when the target function meets the precision requirement, stopping iteration; finally obtaining the anisotropy parameter epsilon of the inversion result reaching the logging resolution(V)、δ(V)、γ；
Wherein,
in the objective function, the extended expression of the Jacobian matrix J in the AVO is as follows:
the first-order partial derivatives of the reflection coefficients with respect to the elastic and anisotropic parameters are expressed by the central difference method:
by applying the joint inversion of PP wave, PS1 wave and PS2 wave AVA, the Jacobian matrix J is as follows:
and H ═ JTJ, then there are:
wherein e isPP、ePS1And ePS2Residual errors of PP wave, PS1 wave, PS2 wave AVO, θ1…θiRepresents each incident angle, R, of the longitudinal wave1…RiRepresents PP wave, PS1 wave and,Reflection coefficient of PS2 wave.
2. The fracture medium PP wave and split PS wave AVO joint inversion method of claim 1,
the specific process of step S1 is:
s1.1, extracting a PP wave AVA gather;
1.1.1) calculating the depth of a reflection point and the corresponding incident angle of different shot points according to the track head information of the seismic channel;
1.1.2) extracting amplitude values of incident angles needing inversion to form a new PP wave AVA gather;
s1.2, extracting an AVA trace set of PS1 waves:
1.2.1) calculating the depth of a reflection point and the corresponding incident angle of different shot points according to the track head information of the seismic channel;
1.2.2) extracting amplitude values of incidence angles needing inversion to form a new PS1 wave AVA gather;
1.2.3) compressing the time of the PS1 wave to the vertical travel time of the PP wave by combining the longitudinal wave velocity to obtain a compressed PS1 wave AVA gather:
tPS1(tPP)＝tPP+u(tPP)；
tpprepresenting the vertical travel time, t, of the PP-waveps1(tpp) Time representing PS1 waves; u (tPP) is the transition moveout;
s1.3, PS2 wave AVA gather extraction
S1.3.1, based onCalculating the speed of the slow transverse wave; wherein γ represents the degree of anisotropic development of the transverse wave velocity, Vs2Indicating the speed, V, of the slow transverse waves1Representing the velocity of the fast transverse wave; then, the depth of a reflection point and the incident angle corresponding to different shot points are obtained according to the information of the PS1 wave in the step S1.2;
s1.3.2, extracting amplitude values of the incident angles to be inverted to form a new PS2 wave AVA gather;
s1.3.3, compressing the time of the PS2 wave to the vertical travel time of the PP wave by combining with the longitudinal wave velocity to obtain a compressed PS2 wave AVA gather:
tps2(tpp)＝tpp+u(tpp)；
tpprepresenting the vertical travel time, t, of the PP-waveps2(tpp) Time representing PS2 waves; u (tPP) is the transition moveout.
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 (5)
| 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 |
| CN115494550B (en) * | 2022-09-19 | 2024-10-11 | 成都晶石石油科技有限公司 | Rapid and slow wave splitting crack prediction method and system based on ray domain transverse wave impedance |
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 |
| US11215721B2 (en) * | 2014-08-19 | 2022-01-04 | Cgg Services Sas | 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 | |
| CN104597490B (en) | Multi-wave AVO reservoir elastic parameter inversion method based on accurate Zoeppritz equations | |
| CN106405651B (en) | Full waveform inversion initial velocity model construction method based on logging matching | |
| WO2017167191A1 (en) | Method and device for processing seismic data | |
| CN111025387B (en) | Pre-stack earthquake multi-parameter inversion method for shale reservoir | |
| CN107894612B (en) | A kind of the sound impedance inversion method and system of Q attenuation by absorption compensation | |
| CN109765616B (en) | Amplitude-preserving wave field continuation correction method and system | |
| CN109188511B (en) | Sand shale thin interbed medium multi-wave AVO joint inversion method | |
| CN111290016B (en) | Full waveform speed modeling inversion method based on geological model constraint | |
| CN110780351B (en) | Longitudinal wave and converted wave prestack joint inversion method and system | |
| CN110531410B (en) | Least square reverse time migration gradient preconditioning method based on direct wave field | |
| CN113031068B (en) | Reflection coefficient accurate base tracking prestack seismic inversion method | |
| CN111175821B (en) | Step-by-step inversion method for anisotropic parameters of VTI medium | |
| CN111025388B (en) | Multi-wave combined prestack waveform inversion method | |
| CN110187389B (en) | AVA inversion method based on thin layer reflection theory | |
| CN109143352B (en) | A kind of anisotropic medium Seismic reflection character establishing equation method | |
| CN111239805B (en) | Block constraint time-lapse seismic difference inversion method and system based on reflectivity method | |
| CN110764146B (en) | Space cross-correlation elastic wave reflection waveform inversion method based on acoustic wave operator | |
| CN109490964B (en) | Improved high-precision AVO elastic parameter fast inversion method | |
| CN111308550B (en) | Anisotropic parameter multi-wave joint inversion method for shale VTI reservoir | |
| CN107918152A (en) | A kind of seismic coherence chromatography imaging method | |
| CN110333534A (en) | A kind of Bayes's time shift AVO inversion method and system based on Biot theory | |
| CN112925022B (en) | Method for predicting anisotropic parameters of shale VTI medium | |
| CN106556862B (en) | PP wave reflection coefficient calculation methods for shale gas exploration AVO technologies | |
| CN111077573A (en) | Method, device and system for determining stratum elastic parameters |
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 |