CN117572495A - Quantitative prediction method for crack scale - Google Patents
Quantitative prediction method for crack scale Download PDFInfo
- Publication number
- CN117572495A CN117572495A CN202311378497.6A CN202311378497A CN117572495A CN 117572495 A CN117572495 A CN 117572495A CN 202311378497 A CN202311378497 A CN 202311378497A CN 117572495 A CN117572495 A CN 117572495A
- Authority
- CN
- China
- Prior art keywords
- wave
- anisotropy
- component
- degree
- transverse wave
- 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 45
- 238000002474 experimental method Methods 0.000 claims abstract description 5
- 239000011159 matrix material Substances 0.000 claims description 28
- 239000013598 vector Substances 0.000 claims description 13
- 239000011229 interlayer Substances 0.000 claims description 9
- 239000011435 rock Substances 0.000 claims description 5
- 238000000926 separation method Methods 0.000 claims description 5
- 238000006073 displacement reaction Methods 0.000 claims description 3
- 238000011161 development Methods 0.000 abstract description 5
- 206010017076 Fracture Diseases 0.000 description 43
- 208000010392 Bone Fractures Diseases 0.000 description 38
- 230000004044 response Effects 0.000 description 12
- 239000011148 porous material Substances 0.000 description 10
- 230000015572 biosynthetic process Effects 0.000 description 7
- 238000004364 calculation method Methods 0.000 description 7
- 238000005755 formation reaction Methods 0.000 description 7
- 230000008569 process Effects 0.000 description 7
- 230000001419 dependent effect Effects 0.000 description 6
- 238000010586 diagram Methods 0.000 description 6
- 230000010287 polarization Effects 0.000 description 6
- 238000004458 analytical method Methods 0.000 description 4
- 239000012530 fluid Substances 0.000 description 4
- 239000010410 layer Substances 0.000 description 4
- 229920006395 saturated elastomer Polymers 0.000 description 4
- 238000005259 measurement Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 239000002245 particle Substances 0.000 description 2
- 239000007787 solid Substances 0.000 description 2
- 238000001228 spectrum Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 2
- 230000004888 barrier function Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 238000004090 dissolution Methods 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000013508 migration Methods 0.000 description 1
- 230000005012 migration Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 239000002356 single layer Substances 0.000 description 1
- 238000010183 spectrum analysis Methods 0.000 description 1
- 238000006467 substitution reaction 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. analysis, for interpretation, for correction
-
- 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. analysis, for interpretation, for correction
- G01V1/30—Analysis
-
- 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. analysis, for interpretation, for correction
- G01V1/30—Analysis
- G01V1/303—Analysis for determining velocity profiles or travel times
-
- 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. analysis, for interpretation, for correction
- 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
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/40—Seismology; Seismic or acoustic prospecting or detecting specially adapted for well-logging
- G01V1/44—Seismology; Seismic or acoustic prospecting or detecting specially adapted for well-logging using generators and receivers in the same well
- G01V1/48—Processing data
- G01V1/50—Analysing data
-
- 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/61—Analysis by combining or comparing a seismic data set with other data
- G01V2210/612—Previously recorded data, e.g. time-lapse or 4D
-
- 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/61—Analysis by combining or comparing a seismic data set with other data
- G01V2210/616—Data from specific type of measurement
- G01V2210/6161—Seismic or acoustic, e.g. land or sea measurements
-
- 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/61—Analysis by combining or comparing a seismic data set with other data
- G01V2210/616—Data from specific type of measurement
- G01V2210/6169—Data from specific type of measurement using well-logging
-
- 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/622—Velocity, density or impedance
- G01V2210/6222—Velocity; travel time
-
- 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/64—Geostructures, e.g. in 3D data cubes
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Remote Sensing (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 quantitative prediction method for crack dimensions, which comprises the following steps: s1: acquiring the degree of transverse wave splitting anisotropy by using logging data and empirical values; the empirical value is measured according to a petrophysical experiment; s2: acquiring a transverse wave splitting anisotropy degree observation value by using the vertical seismic section data; s3: constructing an objective function based on a least square method according to the transverse wave splitting anisotropy degree and the transverse wave splitting anisotropy degree observation value; s4: and solving the objective function, wherein the crack radius corresponding to the optimal solution is the crack scale. The method can successfully and quantitatively predict the crack scale and provide technical support for oil and gas reservoir development.
Description
Technical Field
The invention relates to the technical field of oil and gas reservoir development, in particular to a quantitative prediction method for crack dimensions.
Background
Cracks develop extensively in the crust, resulting from crust stresses resulting from extended geologic structure movements. They can be used as seepage channels for reservoir space and migration, thereby controlling reservoir formation and distribution; at the same time, the cracks may also act as barriers to fluid, impeding fluid seepage. Thus, fractures are of great significance in reservoir prediction. However, the scale of natural fracture development in actual formations tends to vary due to differences in the direction, magnitude, and timing of the earth stresses.
Cracks can be classified into three types according to scale: large scale fractures (over 1/4 wavelength), such as formation scale fractures; mesoscale fractures (between 1/4 wavelength and 1/100 wavelength), such as intra-layer fractures; and microscale cracks (less than 1/100 wavelength), such as dissolution cracks and pore throats of rock grain size. In fractured reservoirs, conditions (such as reservoir pore connectivity) depend not only on fracture direction and fracture density, but also are closely related to fracture dimensions. In particular, for shale or tight sandstone reservoirs, fracture dimensions have an important impact on engineering decisions.
The seismic data contains fracture information and the seismic method is very important for identifying the fracture. When the fracture scale is greater than the seismic resolution (1/4 wavelength), it appears as a fault, bend or fold in the seismic profile. At this time, a series of technical methods such as coherence and curvature methods may be used to identify the crack. However, cracks with dimensions below the seismic resolution (1/4 wavelength) are difficult to identify by such changes in events and must be detected by seismic anisotropy.
There are many crack prediction methods based on the seismic anisotropy theory, and two main types are: prestack azimuth anisotropy analysis methods based on longitudinal wave anisotropy theory, such as amplitude and velocity changes with azimuth, and Shear Wave Splitting (SWS) analysis based on shear wave anisotropy theory. The above method characterizes the relationship between fracture parameters and seismic response based on a fracture equivalent model, and the petrophysical model uses anisotropic parameters, fracture density, and other parameters to describe the extent of fracture development. However, media containing a small number of large-scale cracks or a large number of small-scale cracks may exhibit the same anisotropy. Therefore, the above method cannot quantitatively predict the crack scale.
Disclosure of Invention
The invention aims to provide a quantitative prediction method for crack scale.
The technical scheme of the invention is as follows:
a quantitative prediction method for crack scale comprises the following steps:
s1: acquiring the degree of transverse wave splitting anisotropy by using logging data and empirical values; the empirical value is measured according to a petrophysical experiment;
s2: acquiring a transverse wave splitting anisotropy degree observation value by using the vertical seismic section data;
s3: constructing an objective function based on a least square method according to the transverse wave splitting anisotropy degree and the transverse wave splitting anisotropy degree observation value;
s4: and solving the objective function, wherein the crack radius corresponding to the optimal solution is the crack scale.
Preferably, in step S1, the obtaining the degree of transverse wave splitting anisotropy specifically includes the following substeps:
s11: constructing an elastic stiffness matrix, and acquiring an elastic stiffness tensor according to the elastic stiffness matrix;
s12: substituting the elastic stiffness tensor into a wave propagation equation to obtain an approximate analytic expression of the phase velocity, and calculating to obtain a fast transverse wave velocity and a slow transverse wave velocity according to the approximate analytic expression of the phase velocity;
s13: and calculating and obtaining the anisotropy degree of the transverse wave splitting according to the fast transverse wave speed and the slow transverse wave speed.
Preferably, in step S11, the elastic stiffness matrix is:
wherein: c is the effective elastic tensor; e-shaped article 0 For the applied strain field; c (C) m Is the elasticity of rock matrixTensors; n is the total number of inclusions;volume fraction of the nth inclusion; e-shaped article inc Sum sigma inc The strain field and stress inside each inclusion are respectively; sigma (sigma) 0 Is a stress field;
in step S12, the wave propagation equation is:
wherein: ρ is the medium density; u is a displacement vector; t is a time variable; l is a partial derivative operator matrix; t represents a transpose; f is a physical strength vector;
the approximate analytical formula of the phase velocity is:
A＝-(Γ 11 +Γ 22 +Γ 33 ) (8)
B＝Γ 11 Γ 22 +Γ 11 Γ 33 +Γ 22 Γ 33 -Γ 12 2 -Γ 13 2 -Γ 23 2 (9)
C'＝-(Γ 11 Γ 22 Γ 33 +2Γ 12 Γ 13 Γ 23 -Γ 12 2 Γ 33 -Γ 13 2 Γ 22 -Γ 23 2 Γ 11 ) (10)
wherein: v is the phase velocity; when m=1, v 1 Is the wave velocity of the longitudinal wave; when m=2, v 2 Is the slow transverse wave velocity; when m=3, v 3 Is the fast transverse wave velocity; H. a, B, C ', ψ', Δ, G are all intermediate parameters; Γ -shaped structure 11 、Γ 22 、Γ 33 、Γ 12 、Γ 13 、Γ 23 Tensors in the christiffel matrix, both HTI media;
in step S13, the degree of transverse wave splitting anisotropy is calculated by the following formula:
wherein:is the degree of shear wave splitting anisotropy.
Preferably, tensors in the Christoffel matrix of each HTI medium are respectively:
Γ 23 ＝Γ 32 ＝(c 33 -c 44 )n y n z (15)
wherein: c 11 、c 12 、c 66 、c 13 、c 33 、c 44 Are elastic stiffness tensors in an elastic stiffness matrix of the HTI medium; n is n x 、n y 、n z Are propagation direction vectors n= (n) x ,n y ,n z ) T Is a value of (2).
Preferably, in step S2, the obtaining of the observed value of the degree of transverse wave splitting anisotropy specifically includes the following substeps:
s21: according to the vertical seismic section, obtaining an R component and a T component of a VSP wave;
s22: performing wave field separation according to the R component and the T component to obtain an up-going wave field;
s23: and taking a time window in the range of 100-250ms below the first arrival of the up-going wave field, and calculating and obtaining the observed value of the transverse wave splitting anisotropy degree in the time window.
Preferably, in step S23, the observed value of the degree of transverse wave splitting anisotropy is calculated by the following formula:
wherein:is the observed value of the degree of transverse wave splitting anisotropy; δt is the time delay between the fast and slow shear waves; t is t s Is the interlayer propagation time of the slow shear wave.
Preferably, the time delay between the fast and slow shear waves is calculated by:
P(w)＝RE R (w)RE T (w)+IM R (w)IM T (w) (21)
wherein: f is the frequency;the azimuth angle formed by the crack direction and the R component under the condition of any integer e; p is the real part of the R component multiplied by the real part of the T component, plus the imaginary part of the R component multiplied by the imaginary part of the T component; AM (AM) T Amplitude for the T component; AM (AM) R The amplitude of the R component; p (P) j P is corresponding to the j-th seismic data; q (Q) j Q is corresponding to the j-th seismic data; zeta type toy j An included angle between a measuring line corresponding to the jth seismic data and the radial direction is formed; p (w) is P in the frequency domain; RE (RE) R (w) is the real part of the R component in the frequency domain; RE (RE) T (w) is the real part of the T component in the frequency domain; IM (instant Messaging) R (w) is the imaginary part of the R component in the frequency domain; IM (instant Messaging) T (w) is the imaginary part of the T component in the frequency domain; q (w) is Q in the frequency domain; AM (AM) T (w) is AM in the frequency domain T 。
Preferably, the interlayer propagation time of the slow shear wave is obtained by the sub-steps of:
the fast shear wave and the slow shear wave are separated by:
wherein: s1 (t) is a fast S wave at the moment t; s2 (t- δt) is a slow S wave at the time of t- δt;azimuth formed for the fracture direction and the R component; r (t) is the R component at time t; t (T) is the T component at time T;
the interlayer propagation time of the slow shear wave can be obtained by selecting the propagation time of the separated slow shear wave.
Preferably, in step S3, the objective function is:
wherein: l (alpha) f ) Is an objective function;the degree of transverse wave splitting anisotropy at the ith frequency point;the observed value is the transverse wave splitting anisotropy degree at the ith frequency point; n' is the number of frequency bins.
The beneficial effects of the invention are as follows:
the method establishes quantitative relation between shear wave frequency-related response and different fracture scales when acquiring the degree of shear wave splitting anisotropy, obtains model parameters by using logging data, and obtains theoretical curves of shear wave frequency response of different fracture scales; then, extracting observed data of shear wave frequency response from the actually converted shear wave data; and finally, establishing a least square objective function by utilizing the established quantitative relation and observation data, and obtaining a crack scale by successfully and quantitatively predicting the optimal solution of the least square objective function, thereby providing technical support for the development of the oil and gas reservoir.
Drawings
In order to more clearly illustrate the embodiments of the invention or the technical solutions of the prior art, the drawings which are used in the description of the embodiments or the prior art will be briefly described, it being obvious that the drawings in the description below are only some embodiments of the invention, and that other drawings can be obtained according to these drawings without inventive faculty for a person skilled in the art.
FIG. 1 is a flow chart of a quantitative prediction method of crack dimensions according to the present invention;
FIG. 2 is a schematic diagram of an observation system of a VSP according to an embodiment;
FIG. 3 is a schematic diagram of a VSP data reception process according to an embodiment;
FIG. 4 is a schematic diagram of a process of S-wave generation and acceptance of a synthetic data observation system according to another embodiment;
FIG. 5 is a schematic diagram of a coordinate system at the level of a synthetic data observation system according to another embodiment;
FIG. 6 is a graph showing shear wave anisotropy as a function of frequency for different fracture dimensions for one embodiment;
FIG. 7 is a schematic representation of the R component of synthetic data with an azimuth angle in the range of 0-180 for one embodiment;
FIG. 8 is a schematic diagram of the T component of synthesized data with an azimuth angle in the range of 0-180 for one embodiment;
FIG. 9 is a graph showing observations of shear wave splitting anisotropy in one embodiment;
FIG. 10 is a graph showing the results of objective function values at different fracture scales according to one embodiment;
FIG. 11 is a graph showing shear wave anisotropy as a function of frequency for various fracture dimensions at 4010m depth;
FIG. 12 is a graph showing shear wave anisotropy as a function of frequency for various fracture dimensions at a depth of 4020 m;
FIG. 13 is a schematic diagram of an observation system for actual VSP data in accordance with one embodiment;
FIG. 14 is a graph of the upgoing wave results of the R component of a particular embodiment Y2 well VSP wavefield separation;
FIG. 15 is a graph of the upgoing wave results of the T component of a particular embodiment Y2 well VSP wavefield separation;
FIG. 16 is a graph showing the results of the amplitude spectrum of the R component of one embodiment;
FIG. 17 is a graph showing the amplitude spectrum results of the T component of one embodiment;
FIG. 18 is a graph showing shear wave anisotropy as a function of frequency at various fracture scales for a particular embodiment, and observations of the degree of shear wave splitting anisotropy at different frequency points obtained using VSP data at 4010m depth;
FIG. 19 is a graph showing the results of an objective function at different fracture scales for 4010m depths according to an embodiment;
FIG. 20 is a graph showing the results of crack scale calculations for different depths according to one embodiment.
Detailed Description
The invention will be further described with reference to the drawings and examples. It should be noted that, without conflict, the embodiments and technical features of the embodiments in the present application may be combined with each other. It is noted that all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this application belongs unless otherwise indicated. The use of the terms "comprising" or "includes" and the like in this disclosure is intended to cover a member or article listed after that term and equivalents thereof without precluding other members or articles.
As shown in fig. 1, the invention provides a quantitative prediction method for crack scale, which comprises the following steps:
s1: acquiring the degree of transverse wave splitting anisotropy by using logging data and empirical values; the empirical values are measured according to petrophysical experiments.
In a specific embodiment, the degree of shear wave splitting anisotropy is obtained by:
first, a multi-scale equivalent media model is established. Consider the case of a pore space consisting of a random isotropic collection of microcracks and spherical pores, each crack and pore radius being determined by the grain size. At the same time, a set of aligned cracks with a radius greater than the grain size are introduced, and the microcracks and cracks are assumed to have a common aspect ratio, denoted by r. The volume of individual cracks, pores and fissures is calculated by:
wherein: c v 、p v And f v The volumes of individual cracks, pores and fissures, respectively; a is the radius of the crack and the aperture; a, a f Is the radius of the crack.
In the multi-scale equivalent media model described above, the cracks are larger than the cracks and pores, with each crack and pore connecting at most one crack and no other crack. To achieve this, the number of cracks and voids must far exceed the number of cracks. After these relationships are established, the effective elastic constants of the embedded inclusion materials are calculated using the interaction energy method, specifically by the following formula:
wherein: c is the effective elastic tensor; e-shaped article 0 For the applied strain field; c (C) m Is the elastic tensor of the rock matrix; n is the total number of inclusions;volume fraction of the nth inclusion; e-shaped article inc Sum sigma inc The strain field and stress inside each inclusion are respectively; sigma (sigma) 0 Is a stress field;
the formula (1) is an elastic stiffness matrix, which is obtained by deducting the following formula:
from the pressure field, the explicit equation for the elastic tensor can be expressed as:
wherein: c (C) 11 、C 33 、C 44 、C 12 、C 13 Are elastic tensors in the elastic stiffness matrix; lambda is the first parameter of lame; mu is the second parameter of lame;and->Volume fraction or porosity associated with cracks, pores, and fissures, respectively; l (L) 2 、G 1 、G 2 、G 3 、D 1 、D 2 、F 1 、F 2 、L 4 Are all intermediate parameters; sigma (sigma) c Both kappa are intermediate process variables; v is poisson's ratio;
the specific expression of the intermediate parameters is as follows:
wherein: iota, beta, tau f 、γ、γ' is the intermediate process quantity; i is an imaginary unit; w is the angular frequency; τ m Relaxation time for water saturation of rock; k (K) c Is an intermediate process quantity, and K is smaller due to the small aspect ratio of cracks c ＝0；
The intermediate process quantities are calculated by the following formulas:
in the middle of：∈、∈ f All are intermediate process quantities; gr is the size of the particle size; ρ s And ρ f Densities of saturated solids and fluid, respectively; v (V) s Shear wave velocity for saturated solids; v (V) f Is the speed of sound in the fluid.
According to the formulas (29) - (51), an elastic coefficient matrix (elastic stiffness matrix) of the multi-scale fracture equivalent medium theoretical model shown in the formula (1) can be obtained, and a foundation is provided for further quantifying the anisotropy of the model.
Then, calculation of the degree of SWS anisotropy was performed. When the elastic coefficient matrix of the multiscale fracture equivalent medium theoretical model is known (formula (1)), the phase velocity generated by the seismic waves through the model can be obtained by utilizing the elastic wave theory. The wave equation (wave propagation equation) for any uniform anisotropic medium is as follows:
wherein: ρ is the medium density; u is a displacement vector; t is a time variable; l is a partial derivative operator matrix; t represents a transpose; f is a physical strength vector.
The physical strength ρf is ignored, and the propagation direction vector of the wave is assumed to be n= (n) x ,n y ,n z ) T The position vector is x= (x) x ,x y ,x z ) T The polarization vector is p= (P x ,p y ,p z ) T The method comprises the steps of carrying out a first treatment on the surface of the v is the phase velocity. The required Christoffel equation can be obtained by solving the following equation:
depending on the symmetry of the elastic matrix, the Christiffel matrix is also symmetrical, i.e., Γ 12 ＝Γ 21 ，Γ 13 ＝Γ 31 And Γ 23 ＝Γ 32 . The Christoffel matrix of the horizontal transversely isotropic medium is represented as follows:
Γ 23 ＝Γ 32 ＝(c 33 -c 44 )n y n z (15)
wherein: c 11 、c 12 、c 66 、c 13 、c 33 、c 44 Are elastic stiffness tensors in an elastic stiffness matrix of the HTI medium; n is n x 、n y 、n z Are propagation direction vectors n= (n) x ,n y ,n z ) T Is a value of (2).
The normal unit vector of the incident wavefront is expressed as follows:
wherein: θ is the angle between the propagation direction and the z-axis of the symmetry axis;is the angle between the propagation direction and the x-axis of the symmetry axis.
The Christoffel equation for a horizontal transversely isotropic medium is made non-zero solved by:
by simplification, it can become a one-dimensional cubic equation, so that ρV can be solved 2 ：
ρV 2 ) 3 +A(ρV 2 ) 2 +BρV 2 +C＝0 (55)
Wherein:
A＝-(Γ 11 +Γ 22 +Γ 33 ) (8)
B＝Γ 11 Γ 22 +Γ 11 Γ 33 +Γ 22 Γ 33 -Γ 12 2 -Γ 13 2 -Γ 23 2 (9)
C'＝-(Γ 11 Γ 22 Γ 33 +2Γ 12 Γ 13 Γ 23 -Γ 12 2 Γ 33 -Γ 13 2 Γ 22 -Γ 23 2 Γ 11 ) (10)
substitution into ρV 2 =x, then equation (55) becomes:
x 3 +Ax 2 +Bx+C＝0 (56)
the real root of this third equation of equation (56) is:
wherein:
thus, ρV can be obtained by the formula (57) 2 . When m=1, 2,3, the accurate values of the phase velocities of the qP, qS2, and qS1 waves satisfy qP>qS1>qS2. By means ofTaylor expansion of (2), then->Omitting higher order items, we can get +.>Equation (57) can be approximated as:
the approximate analytical formula of the phase velocity can be obtained as follows:
in the formula (3) -formula (10): v is the phase velocity; when m=1, v 1 Is the wave velocity of the longitudinal wave; when m=2, v 2 Is the slow transverse wave velocity; when m=3, v 3 Is the fast transverse wave velocity; H. a, B, C ', ψ', Δ, G are all intermediate parameters; Γ -shaped structure 11 、Γ 22 、Γ 33 、Γ 12 、Γ 13 、Γ 23 Are tensors in the christiffel matrix of the HTI medium.
The degree of anisotropy of the SWS (degree of transverse wave splitting anisotropy) can be obtained using the slow transverse wave velocity and the fast transverse wave velocity as follows:
wherein:is the degree of shear wave splitting anisotropy.
The degree of SWS anisotropy (degree of shear wave splitting anisotropy) is obtained when the parameters in formulas (29) - (51) are obtained from well log data and empirical values as a function of fracture dimensions and frequency
S2: and obtaining the transverse wave splitting anisotropy degree observation value by using the vertical seismic section data.
When the VSP data is time window selected for the vertical seismic profile data, as shown in fig. 2 and 3. In fig. 3, the source excited at point O produces a longitudinal wave; after reflection from interfaces B and D at depths B and D, converted shear waves are generated and ultimately received by the receiving points at depths a and C, respectively. The (1) wave is the head wave of the longitudinal wave, the (2) wave and the (4) wave are fast shear waves, and the (3) wave and the (5) wave are slow shear waves; f and e are times when the receiving point receives the first arrival of the shear wave, a and c are times when the receiving point receives the uplink fast shear wave, and b and d are times when the receiving point receives the uplink slow shear wave.
Although downstream waves may also be received during the formation of a VSP single record event, downstream main waves do not contain waves reflected from interfaces below the depth of the point of reception. However, the up-going shear waves received at time points a and b may provide information of the structure in the depth segment AB, whereas the down-going shear waves received at time points c and d only contain information of the properties in the depth segment CD. Thus, the VSP data should be split into up-going and down-going waves to obtain the wavefield response of the subsurface medium, and the up-going wavefield should be further analyzed.
Furthermore, to obtain more prominent SWS information, it is necessary to use an appropriate time window for the upwave field, as shown by the first and second lines from left to right in FIG. 2. By processing the seismic information in these two lines, shear wave anisotropy characteristics of the formation at a depth below the corresponding receivers can be obtained. Logging data is necessary to obtain near-well seismic information. Thus, when the time window is open, an appropriate time frame must be selected prior to the first arrival to obtain seismic information near the well. Seismic information near the well may be used to extract SWS anisotropy near the well.
In a specific embodiment, the obtaining of the observed value of the degree of shear wave splitting anisotropy specifically comprises the following sub-steps:
s21: according to the vertical seismic section, obtaining an R component and a T component of a VSP wave;
s22: performing wave field separation according to the R component and the T component to obtain an up-going wave field;
s23: and taking a time window in the range of 100-250ms below the first arrival of the up-going wave field, and calculating and obtaining the observed value of the transverse wave splitting anisotropy degree in the time window.
The receivers of the Vertical Seismic Profile (VSP) are placed closer to the target layer than the receivers of the surface seismic method. Thus, VSP waves propagate through fewer formations and experience less energy loss than surface seismic waves. Thus, VSP waves may have a high signal-to-noise ratio (S/N) and provide extensive information about the target layer. In the above embodiment, the VSP wave is used to estimate the degree of anisotropy of the shear wave, and the obtained observed value of the degree of anisotropy of the transverse wave splitting can make the final result more accurate.
The original polarization of the shear wave before splitting is along the radial (R) component, because the polarization of the SV wave is along the R component. In order to obtain the frequency variation characteristics of SWS from VSP waves, it is necessary to extract SWS parameters in the frequency domain. First, R and transverse (T) component data in the time domain are converted to the frequency domain. Since the original polarization of the shear wave before splitting is along the R component and the polarization of the SV wave is along the R component, parameters P and Q satisfy:
P(w)＝RE R (w)RE T (w)+IM R (w)IM T (w) (21)
wherein the method comprises the steps ofThe calculation is performed by the following formula:
in equation (20), k may take any integer value,there are four solutions in the 0-pi range, and the correct solution can be obtained by calculating the reciprocal of equation (20). At the time of obtaining azimuth +>The time delay between the fast and slow shear waves can then be obtained by the following calculation:
in the formulas (19) - (23): f is the frequency;the azimuth angle formed by the crack direction and the R component under the condition of any integer e; p is the real part of the R component multiplied by the real part of the T component, plus the imaginary part of the R component multiplied by the imaginary part of the T component; AM (AM) T Amplitude for the T component; AM (AM) R The amplitude of the R component; p (P) j P is corresponding to the j-th seismic data; q (Q) j Q is corresponding to the j-th seismic data; zeta type toy j Between the corresponding measuring line of the jth seismic data and the radial directionAn included angle; p (w) is P in the frequency domain; RE (RE) R (w) is the real part of the R component in the frequency domain; RE (RE) T (w) is the real part of the T component in the frequency domain; IM (instant Messaging) R (w) is the imaginary part of the R component in the frequency domain; IM (instant Messaging) T (w) is the imaginary part of the T component in the frequency domain; q (w) is Q in the frequency domain; AM (AM) T (w) is AM in the frequency domain T 。
By using the time delay δt between the fast shear wave and the slow shear wave, an observed value of the shear wave anisotropy degree (observed value of the shear wave splitting anisotropy degree) can be obtained:
wherein: x is the formation thickness through which the shear wave passes;is the observed value of the degree of transverse wave splitting anisotropy; δt is the time delay between the fast and slow shear waves; t is t s Is the interlayer propagation time of the slow shear wave.
Wherein the interlayer propagation time t of the slow shear wave s Obtained by the following sub-steps:
the fast shear wave and the slow shear wave are separated by:
wherein: s1 (t) is a fast S wave at the moment t; s2 (t- δt) is a slow S wave at the time of t- δt;azimuth formed for the fracture direction and the R component; r (t) is the R component at time t; t (T) is the T component at time T;
the interlayer propagation time of the slow shear wave can be obtained by selecting the propagation time of the separated slow shear wave.
S3: and constructing an objective function based on a least square method according to the transverse wave splitting anisotropy degree and the transverse wave splitting anisotropy degree observation value.
In a specific embodiment, the objective function is:
wherein: l (a) f ) Is an objective function;the degree of transverse wave splitting anisotropy at the ith frequency point;the observed value is the transverse wave splitting anisotropy degree at the ith frequency point; n' is the number of frequency bins.
S4: and solving the objective function, wherein the crack radius corresponding to the optimal solution is the crack scale. I.e. when L (a) f ) With optimal solution, a f I.e. the size of the crack dimension.
In a specific embodiment, a synthetic data test is performed, and the measurement geometry of the synthetic data is shown in fig. 4 and 5. After the longitudinal wave is generated at the surface, it propagates downward until it reaches the reflective interface, after which the shear wave generated at the interface is reflected. A horizontal laterally isotropic medium is located above the reflective interface. The S-wave splits into a fast shear wave and a slow shear wave in a horizontal transversely isotropic medium, which are ultimately received by a detector in the well.
Detectors in the well having two horizontal components, an R component and a T component, are used to record the seismic waveform. In the fracture medium, the polarizations of the fast and slow shear waves are along and perpendicular to the fracture plane (f-axis and s-axis in fig. 5), respectively.
To verify the feasibility of the invention, the required parameters are set. The characteristics of the "average" water-saturated sandstone under high (30-40 MPa) effective stress are obtained. The P wave speed is 4090m s -1 S-wave velocity is 2410m s -1 The calculated poisson ratio is 0.24; density of2370kg m -3 The method comprises the steps of carrying out a first treatment on the surface of the Porosity is 16%; lambda is 1.4X10 10 Pa, μ is 2.1X10 10 Pa. The following values were obtained by measuring a synthetic sandstone sample containing cracks of known geometry and orientation: radius of cracks and voids a=2.75x10 -3 m; aspect ratio r=2×10 -4 m. Furthermore, the time scale of water-saturated sandstone is taken as the time scale, i.e., =2×10 -5 s; the sound velocity in water was 1500m/s -1 . Using these measurements as a rough guide, the remaining parameters can be estimated as follows: the particle size was 2X 10 -4 m, crack Density ε of 0.1, crack Density ε f Is 0.05<0.1。
By the above parameters, an intermediate quantity is obtained by calculation according to the formula (34) -formula (51), wherein the calculation result of the partial parameters is γ=14.0, γ' =1.0. And introducing the intermediate quantity into the formulas (29) - (33) to obtain the elastic coefficient matrix of the multi-scale fracture equivalent medium theoretical model. Subsequently, it is assumed that the angle between the propagation direction of the shear wave and the normal direction of the fracture is 70 °, and different fracture dimensions are set. The frequency-dependent response of the degree of SWS anisotropy is obtained by equation (11), and the result is shown in fig. 6.
In addition, a single-layer anisotropic medium model is established, and the fast polarized shear wave is assumed to be irrelevant to frequency; azimuth angle range is 0-180 degrees, equidistant spacing is 5 degrees; the peak frequency of the wave source was 25Hz. The crack scale of the model was set to 1.5m. Using the frequency dependent response of SWS anisotropy (curve corresponding to 1.5m in fig. 6), the time difference between the fast and slow shear waves was found using equation (24). Based on the model, the R and T components in the frequency domain are first synthesized, and then the R and T components in the time domain are obtained using inverse discrete fourier transform. By combining these parameters, the R component and the T component as shown in fig. 7 and 8 are synthesized. Their time interval is 1ms and duration is 1s.
In order to extract the degree of shear wave anisotropy in the R and T components, the R and T components having an azimuth angle of 75 ° are randomly selected, and an observation of the degree of SWS anisotropy (black dots in fig. 9) is extracted using equation (18). Observations obtained using other azimuth angles are also represented by black dots in fig. 9.
In order to predict the size of the fracture scale, the objective function L and the fracture scale a can be obtained by equation (25) using the observed value (black dot in FIG. 9) and the degree of SWS anisotropy (FIG. 6) f The relationship between them is shown in fig. 10. When a is f When=1.5m, the objective function L reaches a minimum, and there is an optimal solution; the crack scale was 1.5m, which is equal to the set crack scale (curve corresponding to 1.5m in fig. 6). This shows that the method of the present invention for predicting crack dimensions is feasible.
In another specific embodiment, an actual data test is performed, and the fracture scale of the actual data is predicted by adopting the quantitative prediction method of the fracture scale according to the actual logging data and the seismic records.
First, the desired geologic parameters were obtained using the log data of the Y2 well of the shale gas zone of the Sichuan basin, as shown in table 1. Other parameter values not shown are empirical values measured in petrophysical experiments, the same as those cited in the previous examples.
Table 1 geological parameter table
By utilizing the parameters, a multi-scale fracture equivalent medium model is constructed at different depths to obtain shear wave frequency related responses under different fracture scales, and an elastic stiffness matrix C in a frequency domain is obtained ijkl (a f ). The response results for the two depths 4010m and 4020m are shown in fig. 11 and 12.
An observation system of actual VSP data is shown in fig. 13. The shot is located at the surface and a detector in the well receives the actual data (VSP). The offset (100 m) can be considered to be approximately zero compared to the depth of the detector (3090-4090 m). Thus, the direction of the converted shear wave is nearly vertical and the distance traveled by the shear wave is approximately equal to the layer thickness, which in turn is equal to the depth difference between the reflection point and the geophone. The distance between two adjacent detectors is 10m and the sampling rate is 1ms.
Subsequently, three-component VSP numbers are usedObservations of the degree of SWS anisotropy were extracted. The up-going wavefield is obtained by separating the up-going and down-going wavefronts of the R and T components using a Y2 well VSP seismic record by a singular value decomposition algorithm, as shown in fig. 14 and 15. The first line from top to bottom in fig. 14 and 15 represents the first arrival. We have selected a time window within 100-250ms of the first arrival, as shown by the second and third lines from top to bottom in fig. 14 and 15, and use the data in this time window to extract observations of the degree of shear wave anisotropy
Next, taking as an example the underground medium information of 4010m depth, the crack scale of the depth point is predicted. In order to obtain the effective frequency band of the seismic information, spectral analysis is performed on the R component and the T component of different depths. The results of this analysis are shown in fig. 16 and 17. In fig. 16 and 17, the improved data means that the reflection coefficient of the VSP is obtained by Gabor deconvolution first, and then the broadband wavelet is convolved with the reflection coefficient to obtain VSP data with higher dominant frequency and wider frequency band. The effective band range of the improved data is about 6.3-39.8Hz (logarithmic frequency 0.8-1.6).
Subsequently, frequency dependent SWS analysis (within 100 Hz) was performed on the data in the time window to obtain a frequency dependent response of the degree of SWS anisotropy, the results being shown as black dots in fig. 18. The crosses in fig. 18 represent the frequency dependent response of the degree of SWS anisotropy obtained using the improved data. In fig. 18, comparison of the original data and the modified data shows that the diversity of the modified data is reduced. Even if the data in the modified data is weaker in the low frequency band, this does not affect the results, since the resolution of the small-scale cracks increases significantly as the dominant frequency increases.
The cross and frequency dependent response relationship in fig. 18 obtained using the least squares method calculation within 100Hz is calculated according to the objective function shown in equation (25). As shown in fig. 19, when the crack scale is 0.19m, the objective function L has an optimal solution; thus, the crack scale for 4010m depth is 0.19m. Also, the fracture scale of other depth points was predicted to obtain the relationship between the fracture scale and the depth point, and the results of the target interval (4010-4090 m) were shown in fig. 20.
In conclusion, the method can quantitatively predict the crack scale. Compared with the prior art, the invention has obvious progress.
The present invention is not limited to the above-mentioned embodiments, but is intended to be limited to the following embodiments, and any modifications, equivalents and modifications can be made to the above-mentioned embodiments without departing from the scope of the invention.
Claims (9)
1. The quantitative prediction method for the crack scale is characterized by comprising the following steps of:
s1: acquiring the degree of transverse wave splitting anisotropy by using logging data and empirical values; the empirical value is measured according to a petrophysical experiment;
s2: acquiring a transverse wave splitting anisotropy degree observation value by using the vertical seismic section data;
s3: constructing an objective function based on a least square method according to the transverse wave splitting anisotropy degree and the transverse wave splitting anisotropy degree observation value;
s4: and solving the objective function, wherein the crack radius corresponding to the optimal solution is the crack scale.
2. The method according to claim 1, wherein in step S1, the step of obtaining the degree of anisotropy of the shear wave splitting comprises the following steps:
s11: constructing an elastic stiffness matrix, and acquiring an elastic stiffness tensor according to the elastic stiffness matrix;
s12: substituting the elastic stiffness tensor into a wave propagation equation to obtain an approximate analytic expression of the phase velocity, and calculating to obtain a fast transverse wave velocity and a slow transverse wave velocity according to the approximate analytic expression of the phase velocity;
s13: and calculating and obtaining the anisotropy degree of the transverse wave splitting according to the fast transverse wave speed and the slow transverse wave speed.
3. The method according to claim 2, wherein in step S11, the elastic stiffness matrix is:
wherein: c is the effective elastic tensor; e-shaped article 0 For the applied strain field; c (C) m Is the elastic tensor of the rock matrix; n is the total number of inclusions;volume fraction of the nth inclusion; e-shaped article inc Sum sigma inc The strain field and stress inside each inclusion are respectively; sigma (sigma) 0 Is a stress field;
in step S12, the wave propagation equation is:
wherein: ρ is the medium density; u is a displacement vector; t is a time variable; l is a partial derivative operator matrix; t represents a transpose; f is a physical strength vector;
the approximate analytical formula of the phase velocity is:
A＝-(Γ 11 +Γ 22 +Γ 33 ) (8)
B＝Γ 11 Γ 22 +Γ 11 Γ 33 +Γ 22 Γ 33 -Γ 12 2 -Γ 13 2 -Γ 23 2 (9)
C′＝(Γ 11 Γ 22 Γ 33 +2Γ 12 Γ 1323 Γ 12 2 Γ 33 -Γ 13 2 Γ 22 -Γ 23 2 Γ 11 ) (10)
wherein: v is the phase velocity; when m=1, v 1 Is the wave velocity of the longitudinal wave; when m=2, v 2 Is the slow transverse wave velocity; when m=3, v 3 Is the fast transverse wave velocity; H. a, B, C ', ψ', Δ, G are all intermediate parameters; Γ -shaped structure 11 、Γ 22 、Γ 33 、Γ 12 、Γ 13 、Γ 23 Tensors in the christiffel matrix, both HTI media;
in step S13, the degree of transverse wave splitting anisotropy is calculated by the following formula:
wherein:is the degree of shear wave splitting anisotropy.
4. The method of claim 3, wherein tensors in Christoffel matrices for each HTI medium are respectively:
Γ 23 ＝Γ 32 ＝(c 33 -c 44 )n y n z (15)
wherein: c 11 、c 12 、c 66 、c 13 、c 33 、c 44 Are elastic stiffness tensors in an elastic stiffness matrix of the HTI medium; n is n x 、n y 、n z Are propagation direction vectors n= (n) x ，n y ，n z ) T Is a value of (2).
5. The method according to claim 1, wherein in step S2, obtaining the observed value of the degree of anisotropy of the shear wave division specifically includes the following sub-steps:
s21: according to the vertical seismic section, obtaining an R component and a T component of a VSP wave;
s22: performing wave field separation according to the R component and the T component to obtain an up-going wave field;
s23: and taking a time window in the range of 100-250ms below the first arrival of the up-going wave field, and calculating and obtaining the observed value of the transverse wave splitting anisotropy degree in the time window.
6. The method according to claim 5, wherein in step S23, the observed value of the degree of shear wave splitting anisotropy is calculated by the following equation:
wherein:is the observed value of the degree of transverse wave splitting anisotropy; δt is the time delay between the fast and slow shear waves; t is t s Is the interlayer propagation time of the slow shear wave.
7. The method of claim 6, wherein the time delay between the fast shear wave and the slow shear wave is calculated by:
P(w)＝RE R (w)RE T (w)+IM R (w)IM T (w) (21)
wherein: f is the frequency;the azimuth angle formed by the crack direction and the R component under the condition of any integer e; p is the real part of the R component multiplied by the real part of the T component, plus the imaginary part of the R component multiplied by the imaginary part of the T component; AM (AM) T Amplitude for the T component; AM (AM) R The amplitude of the R component; p (P) j P is corresponding to the j-th seismic data; q (Q) j Q is corresponding to the j-th seismic data; zeta type toy j An included angle between a measuring line corresponding to the jth seismic data and the radial direction is formed; p (w) is P in the frequency domain; RE (RE) R (w) is the real part of the R component in the frequency domain; RE (RE) T (w) is the real part of the T component in the frequency domain; IM (instant Messaging) R (w) is the imaginary part of the R component in the frequency domain; IM (instant Messaging) T (w) is the imaginary part of the T component in the frequency domain; q (w) is Q in the frequency domain; AM (AM) T (w) is AM in the frequency domain T 。
8. The method of claim 6, wherein the slow shear wave interlayer propagation time is obtained by the sub-steps of:
the fast shear wave and the slow shear wave are separated by:
wherein: s1 (t) is a fast S wave at the moment t; s2 (t- δt) is a slow S wave at the time of t- δt;azimuth formed for the fracture direction and the R component; r (t) is the R component at time t; t (T) is the T component at time T;
the interlayer propagation time of the slow shear wave can be obtained by selecting the propagation time of the separated slow shear wave.
9. The method according to any one of claims 1-8, wherein in step S3, the objective function is:
wherein: l (a) f ) Is an objective function;the degree of transverse wave splitting anisotropy at the ith frequency point; />The observed value is the transverse wave splitting anisotropy degree at the ith frequency point; n' is the number of frequency bins.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311378497.6A CN117572495B (en) | 2023-10-24 | Quantitative prediction method for crack scale |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311378497.6A CN117572495B (en) | 2023-10-24 | Quantitative prediction method for crack scale |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117572495A true CN117572495A (en) | 2024-02-20 |
| CN117572495B CN117572495B (en) | 2024-05-17 |
Family
ID=
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20020159617A1 (en) * | 2001-03-07 | 2002-10-31 | Lin-Ying Hu | Method for gradually deforming an initial object distribution in a heterogeneous medium, generated by simulation of an object type stochastic model, to best adapt it to imposed physical constraints |
| US20140297193A1 (en) * | 2013-03-26 | 2014-10-02 | Cgg Services Sa | Seismic methods and systems employing shallow shear-wave splitting analysis using receiver functions |
| US20180203145A1 (en) * | 2017-01-13 | 2018-07-19 | Cgg Services Sas | Method and apparatus for unambiguously estimating seismic anisotropy parameters |
| US10677947B1 (en) * | 2018-12-10 | 2020-06-09 | China University Of Petroleum (Beijing) | Horizontal fracture prediction method, device and equipment |
| WO2020140803A1 (en) * | 2019-01-02 | 2020-07-09 | 中国石油天然气股份有限公司 | High-angle crack prediction method, computer device, and computer readable storage medium |
| CN115097520A (en) * | 2022-06-09 | 2022-09-23 | 中国人民解放军国防科技大学 | Frequency-dependent anisotropic shear wave splitting algorithm for multiple VSP data |
| CN115903024A (en) * | 2022-12-26 | 2023-04-04 | 成都理工大学 | Shear wave splitting analysis method based on gradient descent method |
| CN116859455A (en) * | 2022-03-28 | 2023-10-10 | 中国石油化工股份有限公司 | Complex reservoir anisotropy parameter prediction method and device, electronic equipment and medium |
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20020159617A1 (en) * | 2001-03-07 | 2002-10-31 | Lin-Ying Hu | Method for gradually deforming an initial object distribution in a heterogeneous medium, generated by simulation of an object type stochastic model, to best adapt it to imposed physical constraints |
| US20140297193A1 (en) * | 2013-03-26 | 2014-10-02 | Cgg Services Sa | Seismic methods and systems employing shallow shear-wave splitting analysis using receiver functions |
| US20180203145A1 (en) * | 2017-01-13 | 2018-07-19 | Cgg Services Sas | Method and apparatus for unambiguously estimating seismic anisotropy parameters |
| US10677947B1 (en) * | 2018-12-10 | 2020-06-09 | China University Of Petroleum (Beijing) | Horizontal fracture prediction method, device and equipment |
| WO2020140803A1 (en) * | 2019-01-02 | 2020-07-09 | 中国石油天然气股份有限公司 | High-angle crack prediction method, computer device, and computer readable storage medium |
| CN116859455A (en) * | 2022-03-28 | 2023-10-10 | 中国石油化工股份有限公司 | Complex reservoir anisotropy parameter prediction method and device, electronic equipment and medium |
| CN115097520A (en) * | 2022-06-09 | 2022-09-23 | 中国人民解放军国防科技大学 | Frequency-dependent anisotropic shear wave splitting algorithm for multiple VSP data |
| CN115903024A (en) * | 2022-12-26 | 2023-04-04 | 成都理工大学 | Shear wave splitting analysis method based on gradient descent method |
Non-Patent Citations (6)
| Title |
|---|
| O. H. AL-HARRASI等: "Fracture characterization using frequency-dependent shear wave anisotropy analysis of microseismic data", GEOPHYSICAL JOURNAL INTERNATIONAL, vol. 185, no. 02, 31 May 2011 (2011-05-31), pages 1059 - 1070 * |
| YANG, YY等: "Fracture prediction based on walkaround 3D3C vertical seismic profiling data: A case study from the Tarim Basin in China", GEOPHYSICS, vol. 87, no. 04, 12 August 2022 (2022-08-12), pages 123 * |
| 张金伟: "基于频变方位各向异性的裂缝型储层反演方法研究", 基础科学辑, no. 01, 31 January 2020 (2020-01-31), pages 011 - 58 * |
| 曹占宁等: "基于岩石物理模型的转换横波分裂在碳酸盐岩裂缝储层预测中的应用", 石油科学通报, vol. 03, no. 01, 15 March 2018 (2018-03-15), pages 11 - 21 * |
| 曹飞: "裂缝性岩石声波参数实验研究及裂缝性储层测井评价", 基础科学辑, no. 09, 30 September 2015 (2015-09-30), pages 011 - 55 * |
| 杨宇勇: "裂缝发育的地震预测方法研究", 基础科学辑, no. 02, 29 February 2020 (2020-02-29), pages 011 - 230 * |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Zhubayev et al. | Ultrasonic velocity and attenuation anisotropy of shales, Whitby, United Kingdom | |
| Wang | Seismic anisotropy in sedimentary rocks, part 2: Laboratory data | |
| Schoenberg et al. | Seismic anisotropy of fractured rock | |
| Leary et al. | Seismic fracture anisotropy in the Earth's crust: An overview | |
| RU2620785C1 (en) | Method of determining the location of the microseismic focus | |
| US20040054478A1 (en) | Method of wave diagnostics of the oil-and-gas-deposit | |
| NO315767B1 (en) | Procedure for Processing Interface-Reflected Cormorants from an Azimuthal-Anisotropic Geological Formation | |
| Barnhoorn et al. | Experimental identification of the transition from elasticity to inelasticity from ultrasonic attenuation analyses | |
| Li et al. | Experimental study of fracture structure effects on acoustic logging data using a synthetic borehole model | |
| Rabbel et al. | Superdeep vertical seismic profiling at the KTB deep drill hole (Germany): seismic close‐up view of a major thrust zone down to 8.5 km depth | |
| Wang et al. | Analysis and estimation of an inclusion-based effective fluid modulus for tight gas-bearing sandstone reservoirs | |
| Clarke et al. | The influence of fluid type on elastic wave velocity and attenuation in volcanic rocks | |
| King et al. | Mapping faults in the laboratory with seismic scattering 1: the laboratory perspective | |
| Vlahovic et al. | Shear-wave splitting and reservoir crack characterization: the Coso geothermal field | |
| CN117572495B (en) | Quantitative prediction method for crack scale | |
| CN117572495A (en) | Quantitative prediction method for crack scale | |
| Dodds et al. | Experimental and theoretical rock physics research with application to reservoirs, seals and fluid processes | |
| Ali et al. | The accuracy of AVA approximations in isotropic media assessed via synthetic numerical experiments: implications for the determination of porosity | |
| Adam et al. | Mudstone P-wave anisotropy measurements with non-contacting lasers under confining pressure | |
| Ma et al. | Frequency-dependent P-wave anelasticity due to multiscale (fractal) heterogeneities in rocks | |
| CN114488302A (en) | In-situ anisotropic ground stress field prediction method and system | |
| Sadidkhouy et al. | Seismic anisotropy in the crust and upper mantle of the Central Alborz Region, Iran | |
| El Sayed et al. | Petrophysical study of Szolnok Formation, Endrod gas field, Hungary | |
| Bhatti et al. | Low-frequency seismic analysis and direct hydrocarbon indicators | |
| Yu et al. | Quantitative prediction of the fracture scale based on frequency-dependent S-wave splitting |
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 |