A method, system, electronic device and storage medium for crack prediction

Abstract

The present invention discloses a fracture prediction method, system, electronic device and storage medium, relating to the technical field of fracture prediction. The method includes: obtaining the travel time information and total travel time of the current formation at each survey line azimuth; obtaining a total travel time calculation model by means of layer-by-layer inversion; based on the total travel time calculation model and the travel time information of the current formation at each survey line azimuth, inversely solving a first objective function regarding the total travel time to determine the first inversion value of the anisotropic parameter; performing isotropic AVO inversion on the elastic parameters in the incident angle gather to obtain the first inversion value of the elastic parameters; based on the first inversion values of the anisotropic parameter and the elastic parameters, using the Gauss-Newton method to inversely solve a second objective function regarding the first inversion value of the elastic parameters and the first inversion value of the anisotropic parameter to obtain the second inversion value of the anisotropic parameter, and predicting fractures. The present invention improves the accuracy of fracture prediction without the need for logging.

Description

Technical Field

[0001] The present invention relates to the technical field of fracture prediction, and particularly to a fracture prediction method, system, electronic device and storage medium. Background Art

[0002] At present, with the increasing demand for oil and gas energy in various countries around the world, unconventional oil and gas reservoirs represented by fractured reservoirs have gradually become new targets for oil and gas exploration and development to replace conventional oil and gas reservoirs. Fractures are information channels for fluids in fractured reservoirs, connecting karst caves and solution pores in the reservoir, and providing important information for reservoir prediction. Therefore, the existence of fractures has an important impact on the distribution of oil and gas.

[0003] Currently, conventional fracture prediction methods mainly include post-stack seismic attributes and related technologies based on the theory of seismic anisotropy. The post-stack seismic attribute analysis technology is relatively mature, such as coherence attribute, curvature attribute, ant tracking, etc., which can be used to characterize large-scale fractures with a length greater than 1 / 4 seismic wavelength. However, affected by the vertical resolution of seismic data, the post-stack attribute technology cannot identify smaller-scale microfractures. The anisotropic characteristics in the underground medium where microfractures develop are significant. Among them, the horizontal transverse isotropy (HTI) medium, which is usually used to describe a group of vertical fractures and has a horizontal symmetry axis, is a typical anisotropic medium. In the HTI medium, the longitudinal wave amplitude, longitudinal wave velocity (travel time), and shear wave splitting characteristics will change regularly with the change of the incident angle and azimuth angle. Therefore, these three types of information can be used to invert the fracture azimuth angle and anisotropic parameters, and then estimate the development direction and relative development density of fractures.

[0004] However, using shear wave splitting characteristics for fracture prediction requires collecting high-quality shear wave data, generally using downhole three-component data or seafloor four-component data. The results obtained by using the longitudinal wave velocity (travel time) for fracture prediction have low vertical resolution. The method using the longitudinal wave amplitude anisotropy characteristics has higher resolution and strong stability, and is a commonly used method at present. Currently, the methods for anisotropic inversion through longitudinal wave amplitude are mainly divided into the method of longitudinal wave amplitude varying with azimuth (Amplitude Versus Offset and Azimuth, AVAZ) and the method of longitudinal wave amplitude varying with offset (Amplitude Versus Offset, AVO). Most inversion methods require the establishment of an initial model. The initial model of anisotropic parameters is usually established based on the isotropic medium hypothesis or rock physics information or logging information. The anisotropic modeling method based on the isotropic medium hypothesis generally uses the isotropic model as the initial model. The accuracy of this model cannot be guaranteed at all, and it is difficult to adapt to the complex conditions of the actual underground medium. The anisotropic modeling method using rock physics parameters and logging data estimates the anisotropic parameters of the wellbore and establishes an initial model of anisotropic parameters along the layer. However, this modeling method can establish an accurate model near the wellbore, but cannot guarantee the model accuracy far from the well. And due to the limitation of the number of wells drilled during the exploration process, especially in the initial stage of exploration, it is difficult to ensure the accuracy of the initial model, thus affecting the results of the final anisotropic parameter inversion.

[0005] For the modeling method of the initial model of anisotropic parameters in the traditional method, the accuracy of the initial model of anisotropic parameters is low. Or to obtain an initial model with higher accuracy, a large amount of logging data is required, the cost is high, and the modeling process is also relatively complex. Otherwise, it will lead to low accuracy of fracture prediction. Summary of the Invention

[0006] The object of the present invention is to provide a fracture prediction method, system, electronic device and storage medium, which improve the accuracy of fracture prediction without logging.

[0007] To achieve the above object, the present invention provides the following solutions:

[0008] A fracture prediction method, the method includes:

[0009] Obtain the travel time information and total travel time of the current formation at each survey line azimuth; the current formation is the oil and gas reservoir to be predicted; the travel time information is the time when the seismic wave travels from the seismic source through the current formation to the geophone; the total travel time is a function of formation parameters, and the formation parameters include: thickness, longitudinal wave velocity, survey line azimuth angle, longitudinal wave anisotropic parameter, transverse-longitudinal wave anisotropic parameter, fracture azimuth angle and longitudinal wave incident angle;

[0010] A total travel time calculation model is obtained by means of layer-by-layer inversion; the total travel time calculation model is a model regarding anisotropic parameters; the anisotropic parameters include a P-wave anisotropic parameter, an SH-P anisotropic parameter, and an SH-wave anisotropic parameter;

[0011] Based on the total travel time calculation model and the travel time information of the current formation at each survey line azimuth, the least squares method is used to inversely solve a first objective function regarding the total travel time, and the optimal anisotropic parameters in the total travel time calculation model are determined as the first inversion values of the anisotropic parameters;

[0012] Isotropic AVO inversion is performed on the elastic parameters in the incident angle gather to obtain the first inversion values of the elastic parameters; the elastic parameters include P-wave velocity, S-wave velocity, and density;

[0013] Based on the first inversion values of the anisotropic parameters and the first inversion values of the elastic parameters, the Gauss-Newton method is used to inversely solve a second objective function regarding the first inversion values of the elastic parameters and the first inversion values of the anisotropic parameters to obtain the second inversion values of the anisotropic parameters;

[0014] The fractures of the current formation are predicted according to the second inversion values of the anisotropic parameters.

[0015] Optionally, the total travel time calculation model is:

[0016]

[0017] Wherein, Is the total travel time of the kth formation at different survey line azimuths, Is the total travel time of the (k - 1)th formation at different survey line azimuths, Is the interlayer travel time of the kth formation at different survey line azimuths;

[0018]

[0019] Wherein, Is the total travel time formula of the (k - 1)th formation, t j,i Is the interlayer travel time of the ith formation at the jth survey line azimuth, i = 1, 2, …, k - 1, j = 1, 2, …, N;

[0020]

[0021]

[0022] d i Is the thickness of the ith formation, θ i Is the P-wave incident angle of the ith formation, And They are all intermediate parameters, is the P-wave velocity of the i-th formation, ε i is the P-wave anisotropy parameter of the i-th formation, δ i is the P-S wave anisotropy parameter of the i-th formation, is the azimuth angle of the j-th survey line, φ i is the fracture azimuth angle of the i-th formation;

[0023]

[0024] where, d k is the thickness of the k-th formation, θ k is the P-wave incident angle of the k-th formation, and are all intermediate parameters under different survey line azimuths,

[0025] is the P-wave velocity of the k-th formation, ε k is the P-wave anisotropy parameter of the k-th formation, δ k is the P-S wave anisotropy parameter of the k-th formation, is the azimuth angle of the j-th survey line, φ k is the fracture azimuth angle of the k-th formation.

[0026] Optionally, the formula of the first objective function is:

[0027]

[0028] where, J(m) is the first objective function, m is the inversion parameter, is the total travel time of the k-th formation at the survey line azimuth of N.

[0029] Optionally, the second objective function is:

[0030]

[0031] where, ΔM is the second objective function, H is the Hessian matrix, μ is a scalar; I is the identity matrix, is the Jacobin matrix, R obs is the observation coefficient matrix, R0 is the initial reflection coefficient matrix;

[0032]

[0033]

[0034] R PP is the P-wave reflection coefficient, E is the elastic parameter matrix, A is the anisotropy parameter matrix, vp0 is the P-wave velocity, v s0 is the S-wave velocity, ρ is the density, ε is the P-wave anisotropy parameter, δ is the P-S wave anisotropy parameter, γ is the S-wave anisotropy parameter, R(θ1) is the reflection coefficient at the incident angle θ1, R(θ2) is the reflection coefficient at the incident angle θ2, R(θ m ) is the reflection coefficient at the incident angle θ m .

[0035] A fracture prediction system, the system comprising:

[0036] A data acquisition module for acquiring travel time information and total travel time at each survey line azimuth of the current formation; the current formation is the oil and gas reservoir to be predicted; the travel time information is the time for seismic waves to reach the geophone from the seismic source through the current formation; the total travel time is a function of formation parameters, the formation parameters including: thickness, P-wave velocity, survey line azimuth angle, P-wave anisotropy parameter, P-S wave anisotropy parameter, fracture azimuth angle, and P-wave incident angle;

[0037] A total travel time calculation model determination module for obtaining a total travel time calculation model by means of layer-by-layer inversion; the total travel time calculation model is a model for anisotropy parameters; the anisotropy parameters include P-wave anisotropy parameter, P-S wave anisotropy parameter, and S-wave anisotropy parameter;

[0038] An anisotropy parameter primary inversion module for inversely solving a first objective function for the total travel time based on the total travel time calculation model and the travel time information at each survey line azimuth of the current formation by using the least squares method, and determining the optimal anisotropy parameters in the total travel time calculation model as the primary inversion values of the anisotropy parameters;

[0039] An elastic parameter primary inversion module for performing isotropic AVO inversion on the elastic parameters in the incident angle gather to obtain the primary inversion values of the elastic parameters; the elastic parameters include P-wave velocity, S-wave velocity, and density;

[0040] An anisotropy parameter secondary inversion module for inversely solving a second objective function for the primary inversion values of the elastic parameters and the primary inversion values of the anisotropy parameters based on the primary inversion values of the anisotropy parameters and the primary inversion values of the elastic parameters by using the Gauss-Newton method to obtain the secondary inversion values of the anisotropy parameters;

[0041] A prediction module for predicting fractures in the current formation according to the secondary inversion values of the anisotropy parameters.

[0042] An electronic device, comprising:

[0043] One or more processors;

[0044] A storage device on which one or more programs are stored;

[0045] When the one or more programs are executed by the one or more processors, the one or more processors are caused to implement the method as described above.

[0046] A storage medium on which a computer program is stored, wherein the computer program, when executed by a processor, implements the method as described above.

[0047] According to the specific embodiments provided by the present invention, the following technical effects are disclosed by the present invention:

[0048] The present invention discloses a crack prediction method, system, electronic device and storage medium. By using a layer-by-layer inversion method, the travel time of the target layer, the anisotropic parameters of the overlying formation and the crack azimuth angle are used to calculate the preliminary anisotropic parameters and crack azimuth angle of the target layer. Then, the obtained longitudinal wave velocity, transverse wave velocity, medium density, and preliminary anisotropic parameters and crack azimuth angle are optimized to obtain the final longitudinal wave velocity, transverse wave velocity, medium density, and anisotropic parameters, so as to further predict the cracks in the target layer. The present invention improves the accuracy of crack prediction without well logging. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for use in the embodiments. Obviously, the drawings in the following description are only some embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts.

[0050] Figure 1 It is a schematic flow chart of the crack prediction method provided in Embodiment 1 of the present invention;

[0051] Figure 2 It is a schematic diagram of the ray tracing trajectory of a multi-layer anisotropic medium of the present invention;

[0052] Figure 3 It is a system framework diagram of the travel time anisotropy modeling and inversion method based on HTI media;

[0053] Figure 4 It is a schematic diagram of a horizontal layered model with 5 formations;

[0054] Figure 5 It is a schematic diagram of a pre-stack azimuth gather;

[0055] Figure 6 It is a schematic diagram of the in-phase axis at the bottom of the first layer;

[0056] Figure 7Schematic diagram of the in-phase axis at the bottom of the second layer;

[0057] Figure 8 Schematic diagram of the in-phase axis at the bottom of the third layer;

[0058] Figure 9 Schematic diagram of the in-phase axis at the bottom of the fourth layer;

[0059] Figure 10 Schematic diagram of the incident angle gather with azimuth angles of 0°, 60°, and 120°;

[0060] Figure 11 Schematic diagram of the azimuth traveltime anisotropic inversion modeling result;

[0061] Figure 12 Schematic diagram of the isotropic inversion result;

[0062] Figure 13 Schematic diagram of the anisotropic inversion result Figure 1 ;

[0063] Figure 14 Schematic diagram of the anisotropic inversion result Figure 2 ;

[0064] Figure 15 Schematic diagram of the structure of the fracture prediction system provided in Embodiment 2 of the present invention. Detailed implementation manners

[0065] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0066] The purpose of the present invention is to provide a fracture prediction method, system, electronic device, and storage medium, aiming to improve the accuracy of fracture prediction without well logging, and can be applied to the field of fracture prediction technology.

[0067] To make the above objects, features, and advantages of the present invention more obvious and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific implementation manners.

[0068] Embodiment 1

[0069] Figure 1 Schematic diagram of the flow of the fracture prediction method provided in Embodiment 1 of the present invention. As Figure 1 shown, the fracture prediction method in this embodiment includes:

[0070] Step 101: Obtain the travel time information and total travel time of the current formation at each survey line azimuth; the current formation is the oil and gas reservoir to be predicted; the travel time information is the time for the seismic wave to travel from the seismic source through the current formation to the geophone; the total travel time is a function of formation parameters, and the formation parameters include: thickness, P-wave velocity, survey line azimuth angle, P-wave anisotropy parameter, P-S wave anisotropy parameter, fracture azimuth angle, and P-wave incident angle.

[0071] Step 102: Obtain the total travel time calculation model by means of layer-by-layer inversion; the total travel time calculation model is a model about anisotropy parameters; the anisotropy parameters include P-wave anisotropy parameter, P-S wave anisotropy parameter, and S-wave anisotropy parameter.

[0072] Step 103: Based on the total travel time calculation model and the travel time information of the current formation at each survey line azimuth, use the least squares method to inversely solve the first objective function regarding the total travel time, and determine the optimal anisotropy parameters in the total travel time calculation model as the primary inversion values of the anisotropy parameters.

[0073] Step 104: Perform isotropic AVO inversion on the elastic parameters in the incident angle gather to obtain the primary inversion values of the elastic parameters; the elastic parameters include P-wave velocity, S-wave velocity, and density.

[0074] Step 105: Based on the primary inversion values of the anisotropy parameters and the primary inversion values of the elastic parameters, use the Gauss-Newton method to inversely solve the second objective function regarding the primary inversion values of the elastic parameters and the primary inversion values of the anisotropy parameters to obtain the secondary inversion values of the anisotropy parameters.

[0075] Step 106: Predict the fractures of the current formation according to the secondary inversion values of the anisotropy parameters.

[0076] As an optional implementation manner, the total travel time calculation model is:

[0077]

[0078] where is the total travel time of the k-th formation at different survey line azimuths, is the total travel time of the (k - 1)-th formation at different survey line azimuths, is the interlayer travel time of the k-th formation at different survey line azimuths.

[0079]

[0080] where is the total travel time formula of the (k - 1)-th formation, t j,i is the interlayer travel time of the i-th formation at the j-th survey line azimuth, i = 1, 2, …, k - 1, j = 1, 2, …, N.

[0081]

[0082]

[0083] d i is the thickness of the i-th formation, θ i is the incident angle of the P-wave of the i-th formation, and are both intermediate parameters, is the P-wave velocity of the i-th formation, ε i is the P-wave anisotropy parameter of the i-th formation, δ i is the P-S wave anisotropy parameter of the i-th formation, is the azimuth of the j-th survey line, φ i is the fracture azimuth of the i-th formation.

[0084]

[0085] wherein, d k is the thickness of the k-th formation, θ k is the incident angle of the P-wave of the k-th formation, and are both intermediate parameters under different survey line azimuths,

[0086] is the P-wave velocity of the k-th formation, ε k is the P-wave anisotropy parameter of the k-th formation, δ k is the P-S wave anisotropy parameter of the k-th formation, is the azimuth of the j-th survey line, φ k is the fracture azimuth of the k-th formation.

[0087] As an optional implementation manner, the formula of the first objective function is:

[0088]

[0089] wherein, J(m) is the first objective function, m is the inversion parameter, is the total travel time of the k-th formation at the survey line azimuth of N.

[0090] As an optional implementation manner, the second objective function is:

[0091]

[0092] wherein, ΔM is the second objective function, H is the Hessian matrix, μ is a scalar; I is the identity matrix, is the Jacobin matrix, and R obs is the observation coefficient matrix, and R0 is the initial reflection coefficient matrix.

[0093]

[0094]

[0095] R PP is the P-wave reflection coefficient, E is the elastic parameter matrix, A is the anisotropic parameter matrix, v p0 is the P-wave velocity, v s0 is the S-wave velocity, ρ is the density, ε is the P-wave anisotropic parameter, δ is the P-S wave anisotropic parameter, γ is the S-wave anisotropic parameter, R(θ1) is the reflection coefficient at the incident angle θ1, R(θ2) is the reflection coefficient at the incident angle θ2, and R(θ m ) is the reflection coefficient at the incident angle θ m .

[0096] Specifically, the fracture prediction method based on travel-time anisotropy modeling specifically includes:

[0097] I. Azimuth travel-time anisotropy modeling

[0098] (1) Seismic wave travel-time picking

[0099] Using the travel-time information in the pre-stack azimuth angle gather for azimuth travel-time anisotropy inversion, first, the travel-time information of the seismic waves in the target layer needs to be obtained (the travel-time information of the target layer refers to the time from the seismic source to the target layer and then from the target layer to the geophone, that is, twice the time passed from the seismic source to the target layer; the target layer generally refers to the oil and gas-bearing reservoir, which needs to be determined according to specific exploration tasks. Fractures can store oil and gas or provide channels for the migration of oil and gas). In the selected time window, a travel-time scanning function is used to obtain the seismic wave travel-time:

[0100]

[0101] In the formula, is the azimuth angle of the survey line, and N = 1, 2, 3,... is the survey line number; represents the amplitude value at time T within the time window; T0 and T1 respectively represent the lower and upper limits of the selected time window.

[0102] Within the selected time window, when the amplitude value is the maximum, the travel-time scanned at this time is the travel-time of the target layer in the azimuth of this survey line.

[0103] (2) Travel-time of layered anisotropic media

[0104] According to Figure 2The ray-tracing trajectory of the multi-layered anisotropic medium shown (i.e., a process diagram of seismic waves from the source to the target layer and then received by the geophone. According to Figure 2 some geometric relationships of the seismic wave ray trajectory can be known, so as to obtain the travel distance of the ray trajectory), the P-wave travel time in the HTI medium under azimuthal anisotropy is:

[0105]

[0106] Where:

[0107] In the formula: k = 1, 2,..., n is the formation serial number (here the formation refers to several layers from the surface to the underground, and the target layer is only one or several of them); is the P-wave velocity at vertical incidence in the k-th layer; ε k , δ k are the anisotropic parameters of the k-th layer (there are three anisotropic parameters in total, all of which can represent the degree of fracture development. Among them, ε is an anisotropic parameter related to the P-wave; δ is a transitional anisotropic parameter related to the P-wave and S-wave; γ is an anisotropic parameter related to the S-wave. The subscript represents the anisotropic parameter of the k-th layer); θ k is the P-wave incident angle in the k-th layer; d k is the formation thickness of the k-th layer; φ k is the azimuth angle of the fracture symmetry axis of the k-th layer, and are intermediate variables, which are for convenient calculation and have no practical meaning. The P-wave velocity Vp0 and the depth d are known quantities, and the others are unknown quantities.

[0108] For the P-wave velocity it can be calculated according to the travel-time relationship using the Dix formula, that is:

[0109]

[0110] In the formula, T 0,k is the zero-offset time from the first layer to the k-th layer, T 0,k-1 is the zero-offset time from the first layer to the k - 1-th layer, v R,k is the root-mean-square velocity from the first layer to the k-th layer, all of which can be obtained from the seismic record.

[0111] θ1, θ2,..., θ k determine the seismic wave propagation path and satisfy Fermat's principle. Therefore, its constraint conditions are:

[0112]

[0113]

[0114] where X is the offset; λ is the Lagrangian operator, d i is the thickness of the i-th formation, and θ i is the P-wave incident angle of the i-th formation.

[0115] Equation (2) shows that the travel time of the target layer is affected by the anisotropy of the overlying formations and the fracture azimuth. Therefore, in order to calculate the anisotropic parameters of the target layer, the anisotropic parameters of the overlying formations and the fracture azimuth should be obtained first.

[0116] Therefore, an iterative inversion method layer by layer is adopted to calculate the anisotropic parameters and fracture azimuth of the target layer by using the travel time of the target layer, the anisotropic parameters of the overlying formations and the fracture azimuth.

[0117] As Figure 3 shown, first, the parameters of the first layer are obtained by using the travel time in the single-layer case (n = 1 in Equation (2)). Then, the parameters of the underlying formations are obtained layer by layer.

[0118] According to Equation (2), the interlayer travel time of the k-th layer is obtained as:

[0119]

[0120] where d k is the thickness of the k-th formation, and θ k is the P-wave incident angle of the k-th formation.

[0121] When obtaining the parameters of the k-th layer, the anisotropic parameters and fracture azimuth of the overlying formations of the k-th layer are used as known quantities, and an iterative inversion method layer by layer is adopted, that is, first invert the first layer, then the second layer, and then the third layer. When inverting the k-th layer, the anisotropic parameters of the formations above the k-th layer have been inverted.

[0122] Further derived from Equation (6), the seismic wave travel time of the k-1 layer at different survey line azimuths can be expressed as:

[0123]

[0124] where is the total travel time of the k-1 formation, and t j,i is the interlayer travel time of the i-th formation at the j-th survey line azimuth, i = 1, 2,..., k - 1, j = 1, 2,..., N.

[0125]

[0126]

[0127] di is the thickness of the i-th formation, θ i is the incident angle of the P-wave in the i-th formation, and are both intermediate parameters, is the P-wave velocity of the i-th formation, ε i is the P-wave anisotropy parameter of the i-th formation, δ i is the P-S wave anisotropy parameter of the i-th formation, is the azimuth of the j-th survey line, φ i is the fracture azimuth of the i-th formation. When k = 1, t j,i = 0. In the present invention, there are two time concepts. One is t, which is the interlayer travel time between formations, and the other is T, which is the total time for seismic waves to travel from the source to the target formation and then to the geophone. Here, it refers to the interlayer travel time of each formation under different survey line azimuths.

[0128] θ in formula (6) k is unknown before the anisotropy parameters are inverted, that is, Figure 2 the trajectory shown is uncertain. Therefore, the purpose of formula (9) is to constrain Figure 2 the ray trajectory shown, that is, under a specific θ k , it satisfies the principle of the shortest travel time. While iterating the anisotropy parameters in formula 11, θ k is also iterated. Formula (9) is actually a simpler representation of formula (4) and formula (5).

[0129] Since the anisotropy parameters and fracture azimuth of the k-th layer are unknown, the interlayer travel time of seismic waves in the k-th layer under different survey line azimuths is:

[0130]

[0131] Under the constraints of formula (4) and formula (5), the total travel time of seismic waves from the shot point to the geophone under different survey line azimuths is:

[0132]

[0133] According to formula (1), the travel time information of the k-th layer picked up from the prestack azimuth gather under different survey line azimuths is:

[0134]

[0135] In the formula, the superscript T is the transpose, is the total travel time of the k-th layer when the survey line azimuth is , and x = 1, 2…, N.

[0136] According to formula (9) and formula (10), ε can be inverted using the least squares methodk and δ k and φ k . Its objective function (which is the least - squares method, that is, by iteratively calculating the anisotropic parameters until the calculated result is consistent with the actually observed result or the error is minimized) is:

[0137]

[0138] where m is the inversion parameter.

[0139] The parameter γ cannot be directly obtained by the inversion method of P - wave azimuth travel - time anisotropy. However, generally, the variation trends of ε and γ are consistent. Therefore, it can be assumed that there is a linear relationship between them (this relationship can be obtained by linear fitting based on well - logging data).

[0140] γ = Aε + B (12)

[0141] where both A and B are constants.

[0142] By using the least - squares method to iteratively calculate the anisotropic parameters ε k and δ k , the fracture azimuth angle φ k and the seismic - wave incident angle θ of each formation k until the calculated result satisfies formula (11), and the iteration ends. At this time, the anisotropic parameters, fracture azimuth angle, and seismic - wave incident angle of each formation that meet the requirements are obtained. When the synthetic travel - time (the synthetic travel - time is the one with the minimum error according to formula (9) and the actual travel - time , output the anisotropic parameters ε k and δ k and the fracture azimuth angle φ k . Finally, calculate γ according to the linear relationship of formula (12).

[0143] Each imaging point can independently invert the anisotropic parameters using travel - time information and does not require the constraints of horizon and well - logging information. However, since the formation velocity calculated according to the travel - time relationship is an approximate solution, using the layer - by - layer inversion method will inevitably result in the accuracy of the inverted parameters not meeting the requirements of reservoir prediction. To further improve the accuracy of the inverted parameters, the formation anisotropic parameter results obtained by travel - time inversion are used to provide an initial model for subsequent AVAZ inversion.

[0144] II. Anisotropic Inversion

[0145] Based on the assumption of weak anisotropy theory, Rüger derived the formula for the variation of the reflection coefficient R of the longitudinal wave with the incident angle and azimuth angle in HTI media:

[0146]

[0147]

[0148] In the formula, v p0 is the longitudinal wave velocity; v s0 is the shear wave velocity; ρ is the medium density; Z = ρv p0 is the vertical longitudinal wave impedance; is the vertical shear modulus of the shear wave; the superscript "-" represents the average value of the physical quantities of the media above and below the reflection interface; ε, δ, and γ are anisotropy parameters; "Δ" represents the difference in physical quantities of the media above and below the reflection interface; θ is the incident angle of the longitudinal wave; φ is the azimuth angle of the crack symmetry axis; is the azimuth angle of the survey line. (Among them, the azimuth angle φ of the crack symmetry axis, the azimuth angle of the survey line, and the incident angle θ are all known quantities, and the rest are unknown quantities and iterative calculations are required)

[0149] In anisotropic inversion, the Gauss-Newton method is used for anisotropic inversion. In order to improve the inversion accuracy and solve the local minimum solution problem in the Gauss-Newton method, the anisotropic inversion is divided into two stages. At the same time, considering the characteristic that the anisotropy becomes more obvious with the increase of the incident angle of the reflection coefficient. Therefore, in the first stage, an initial elastic parameter model is established based on the elastic parameter data in the well, and the isotropic AVO inversion is performed using the incident angle gather with a small angle (incident angle < 20°) to obtain the elastic parameters (v p0 、v s0and ρ) initial models (the core equation for inversion is Equation (13), and the iterative algorithm is Equation (14). As long as these two equations are known, programming can be carried out for inversion); in the second stage, under the constraints of the initial elastic parameter model and the initial anisotropic parameter model, the anisotropic AVAZ inversion is performed using the mid- and large-angle (incident angle > 20°) incident angle gathers. The purpose of the first stage: Before starting the AVO inversion, only the elastic parameter curve wellhead is known, and the elastic parameters of the entire work area are not known. This stage is to establish an initial elastic parameter model with not very high accuracy based on the elastic parameters wellhead, and obtain a more accurate elastic parameter model for the entire work area through isotropic inversion. The second stage is to perform anisotropic inversion under the constraints of a relatively accurate initial anisotropic parameter model and a relatively accurate initial elastic parameter model, and finally, more accurate elastic parameter and anisotropic parameter results can be obtained. However, this process is often cumbersome and time-consuming. Therefore, in order to optimize the iterative efficiency, the damped least squares method is used to regularize the Gauss-Newton inversion algorithm, and the Euclidean distance similarity is used to establish the stopping condition for iterative inversion, which greatly improves the calculation efficiency and obtains more accurate inversion results. The objective function of the optimization method can be expressed as:

[0150]

[0151] where

[0152]

[0153]

[0154] In the formula, M is the parameter model matrix; μ is a scalar; I is the identity matrix; H is the Hessian matrix; is the Jacobin matrix; R obs and R0 are the observation coefficient matrix and the initial reflection coefficient matrix respectively; R PP is the P-wave reflection coefficient matrix, which can be obtained according to Equation 13; E is the elastic parameter matrix; A is the anisotropic parameter matrix. R(θ1), R(θ2), ……, R(θ m ) are the P-wave reflection coefficients at different P-wave incident angles, which can be obtained according to Equation 13. Robs is known, R0 needs to be iterated each time and is constantly updated, and the remaining quantities are constantly changing during the iteration process

[0155] During the iterative calculation of the elastic parameters (v p0 , v s0 and ρ) and the three anisotropic parameters (ε, δ, and γ), under the given iterative constraints, that is, when R obs is close to R0 (R0 is the initial reflection coefficient matrix, and each time it is iterated, this value is updated using Robs Update once until R0 is consistent with Robs or the error is minimized and the iteration stops), the algorithm iteration stops. Finally, the accurate v can be obtained p0 and v s0 and ρ, ε, δ and γ results. Usually, where fractures develop, the degree of anisotropy is high. Therefore, the fractures in the formation can be predicted based on the finally obtained anisotropic parameter results.

[0156] The following combines specific embodiments to verify the fracture prediction method of the present invention.

[0157] I. Model test

[0158] As Figure 4 shown, a horizontal layered model with 5 formations is established. Among them, the 1st and 5th layers are isotropic layers, and the 2nd, 3rd, and 4th layers are anisotropic layers. Tsuneyama et al. obtained A = 1.2006 and B = -0.0282 in Equation (12) by analyzing the logging data of brine-saturated sandstone and shale. Set γ according to the relationship shown in Equation (13) and also use this relationship in subsequent actual data tests. According to Thomsen's research on anisotropy, most sedimentary rocks are weakly anisotropic, and the value range of their anisotropic parameters is generally 0 to 0.2. The elastic parameters and anisotropic parameters of the model are shown in Table 1.

[0159] Table 1 Physical property parameters and anisotropic parameters of the theoretical model

[0160]

[0161]

[0162] For the prestack azimuth angle gather used for azimuth traveltime anisotropy inversion, set the shot-receiver distance to 280m and the interval of the survey line azimuth angle to 15°; for the incident angle gather used for anisotropy inversion, set the survey line azimuths to 0°, 60°, and 120° and the incident angles to 5° to 40°. According to the above relationship, calculate the PP-wave reflection coefficient of the model shown in Figure 4 in the PP time domain. Then, convolve the reflection coefficient sequence with a Ricker wavelet with a main frequency of 30Hz to generate the prestack azimuth angle gather and incident angle gather shown in Figures 5 - 10 .

[0163] Extract the traveltime of the event corresponding to the bottom interface of each formation. To verify the effectiveness of traveltime inversion, perform a least-squares scan on all possible values of the anisotropic parameters of the first, second, and third layers. The inversion results are shown in Tables 2 - 4. The inversion results show that the maximum error of traveltime inversion is within 10%. For establishing the initial model of anisotropic parameters, this inversion accuracy meets the requirements.

[0164] Analysis of the inversion error results of the first-layer theoretical model in Table 2

[0165]

[0166] Analysis of the inversion error results of the second-layer theoretical model in Table 3

[0167]

[0168] Analysis of the inversion error results of the third-layer theoretical model in Table 4

[0169]

[0170]

[0171] The model test results show that: without any constraint conditions, the travel-time inversion can still obtain relatively accurate results, indicating that the travel-time inversion has good stability.

[0172] Using Equation (12), calculate γ according to the inverted ε and use the final inversion results (ε, δ, γ) to establish the initial model of anisotropic parameters for anisotropic inversion, as Figure 11 shown.

[0173] The anisotropic inversion is divided into two stages. The first stage is isotropic inversion, and the second stage is anisotropic inversion. The isotropic AVO inversion is performed on the multi-layer HTI medium model, and a total of 20 iterations are carried out. When the iteration stops, the residual of the isotropic inversion is very small, about 0.16. The results obtained by the isotropic inversion match well with the true model, as Figure 12 shown.

[0174] Under the constraints of the initial model of anisotropic parameters and the initial model of elastic parameters established by the azimuth travel-time anisotropic inversion method and the isotropic inversion method, the anisotropic inversion is performed on this model, and a total of 25 iterations are carried out. When the iteration stops, the residual of the anisotropic inversion is very small, about 0.32. Comparing the inversion values of the anisotropic parameters with the theoretical values, it can be seen that when there are small errors in the initial model, the inversion values always coincide with the theoretical values, as Figure 13 and Figure 14 shown. Therefore, when the azimuth travel-time anisotropic inversion can provide a relatively accurate initial model of anisotropic parameters, the anisotropic AVO inversion based on the HTI medium can obtain more accurate Thomsen parameters.

[0175] Example 2

[0176] Figure 15 It is a schematic structural diagram of the fracture prediction system provided by Example 2 of the present invention. As Figure 15As shown in the figure, the fracture prediction system in this embodiment includes:

[0177] A data acquisition module 201, configured to acquire the travel time information and the total travel time of the current formation at each survey line azimuth; the current formation is the oil and gas reservoir to be predicted; the travel time information is the time when the seismic wave travels from the seismic source through the current formation to the geophone; the total travel time is a function of formation parameters, and the formation parameters include: thickness, longitudinal wave velocity, survey line azimuth angle, longitudinal wave anisotropy parameter, transverse and longitudinal wave anisotropy parameter, fracture azimuth angle, and longitudinal wave incident angle.

[0178] A total travel time calculation model determination module 202, configured to obtain a total travel time calculation model by means of layer-by-layer inversion; the total travel time calculation model is a model about anisotropy parameters; the anisotropy parameters include longitudinal wave anisotropy parameter, transverse and longitudinal wave anisotropy parameter, and shear wave anisotropy parameter.

[0179] An anisotropy parameter primary inversion module 203, configured to inversely solve a first objective function about the total travel time based on the total travel time calculation model and the travel time information of the current formation at each survey line azimuth by using the least squares method, and determine the optimal anisotropy parameter in the total travel time calculation model as the primary inversion value of the anisotropy parameter.

[0180] An elastic parameter primary inversion module 204, configured to perform isotropic AVO inversion on the elastic parameters in the incident angle gather to obtain the primary inversion value of the elastic parameters; the elastic parameters include longitudinal wave velocity, shear wave velocity, and density.

[0181] An anisotropy parameter secondary inversion module 205, configured to inversely solve a second objective function about the primary inversion value of the elastic parameters and the primary inversion value of the anisotropy parameter based on the primary inversion value of the anisotropy parameter and the primary inversion value of the elastic parameters by using the Gauss-Newton method, and obtain the secondary inversion value of the anisotropy parameter.

[0182] A prediction module 206, configured to predict the fractures of the current formation according to the secondary inversion value of the anisotropy parameter.

[0183] Embodiment 3

[0184] An electronic device includes:

[0185] One or more processors.

[0186] A storage device, on which one or more programs are stored.

[0187] When the one or more programs are executed by the one or more processors, the one or more processors are caused to implement the method as in Embodiment 1.

[0188] Embodiment 4

[0189] A storage medium stores a computer program thereon, wherein the computer program, when executed by a processor, implements the method in Embodiment 1.

[0190] The various embodiments in this specification are described in a progressive manner. Each embodiment focuses on the differences from other embodiments. For the same or similar parts among the various embodiments, reference can be made to each other. For the system disclosed in the embodiments, since it corresponds to the method disclosed in the embodiments, the description is relatively simple. For the relevant parts, reference can be made to the description in the method section.

[0191] Specific examples are used in this article to elaborate on the principles and implementation manners of the present invention. The description of the above embodiments is only used to help understand the method of the present invention and its core idea. At the same time, for those of ordinary skill in the art, according to the idea of the present invention, there will be changes in the specific implementation manners and application scopes. In summary, the content of this specification should not be construed as a limitation on the present invention.

Claims

1. A method for crack prediction, characterized in that, The method includes: Obtaining the travel time information and the total travel time of the current formation at each survey line azimuth; the current formation is the oil and gas reservoir to be predicted; the travel time information is the time for seismic waves to travel from the seismic source through the current formation to the geophone; the total travel time is a function of formation parameters, and the formation parameters include: thickness, P-wave velocity, survey line azimuth angle, P-wave anisotropy parameter, P-S wave anisotropy parameter, fracture azimuth angle, and P-wave incident angle; Obtaining the total travel time calculation model by means of layer-by-layer inversion; the total travel time calculation model is a model of anisotropy parameters; the anisotropy parameters include P-wave anisotropy parameter, P-S wave anisotropy parameter, and S-wave anisotropy parameter; Based on the total travel time calculation model and the travel time information of the current formation at each survey line azimuth, using the least squares method to inversely solve the first objective function regarding the total travel time, and determining the optimal anisotropy parameters in the total travel time calculation model as the first inversion values of the anisotropy parameters; Performing isotropic AVO inversion on the elastic parameters in the incident angle gather to obtain the first inversion values of the elastic parameters; the elastic parameters include P-wave velocity, S-wave velocity, and density; Based on the first inversion values of the anisotropy parameters and the first inversion values of the elastic parameters, using the Gauss-Newton method to inversely solve the second objective function regarding the first inversion values of the elastic parameters and the first inversion values of the anisotropy parameters to obtain the second inversion values of the anisotropy parameters; Predicting the fractures of the current formation according to the second inversion values of the anisotropy parameters.

2. The crack prediction method according to claim 1, characterized in that The total travel time calculation model is: Among them, is the total travel time of the k-th formation at different survey line azimuths, is the total travel time of the (k - 1)-th formation at different survey line azimuths, is the interlayer travel time of the k-th formation at different survey line azimuths; Among them, is the total travel time formula for the (k - 1)-th formation, t j,i is the interlayer travel time of the i-th formation in the j-th survey line azimuth, where i = 1, 2, …, k - 1 and j = 1, 2, …, N; d i is the thickness of the i-th formation, θ i is the incident angle of the longitudinal wave in the i-th formation, and are both intermediate parameters, is the longitudinal wave velocity of the i-th formation, ε i is the longitudinal wave anisotropy parameter of the i-th formation, δ i is the transverse-longitudinal wave anisotropy parameter of the i-th formation, is the azimuth angle of the j-th survey line, φ i is the fracture azimuth angle of the i-th formation; where d k is the thickness of the k-th formation, θ k is the incident angle of the longitudinal wave of the k-th formation, and are both intermediate parameters under different survey line azimuths, is the P-wave velocity of the k-th formation, ε k is the P-wave anisotropy parameter of the k-th formation, δ k is the S-wave to P-wave anisotropy parameter of the k-th formation, is the azimuth of the j-th survey line, φ k is the fracture azimuth of the k-th formation.

3. The crack prediction method according to claim 2, wherein The formula of the first objective function is: Among them, J(m) is the first objective function, and m is the inversion parameter. is the total travel time of the k-th formation at the survey line azimuth of N.

4. The crack prediction method according to claim 1, wherein The second objective function is: where ΔM is the second objective function, H is the Hessian matrix, μ is a scalar; I is the identity matrix, is the Jacobin matrix, R obs is the observation coefficient matrix, and R0 is the initial reflection coefficient matrix; R PP is the longitudinal wave reflection coefficient, E is the elastic parameter matrix, A is the anisotropic parameter matrix, v p0 is the longitudinal wave velocity, v s0 is the shear wave velocity, ρ is the density, ε is the longitudinal wave anisotropic parameter, δ is the shear-longitudinal wave anisotropic parameter, γ is the shear wave anisotropic parameter, R(θ1) is the reflection coefficient at the incident angle θ1, R(θ2) is the reflection coefficient at the incident angle θ2, R(θ m ) is the reflection coefficient at the incident angle θ m .

5. A crack prediction system, characterized in that, The system includes: A data acquisition module, configured to obtain the travel time information and the total travel time of the current formation at each survey line azimuth; the current formation is the oil and gas reservoir to be predicted; the travel time information is the time for seismic waves to travel from the seismic source through the current formation to the geophone; the total travel time is a function of formation parameters, and the formation parameters include: thickness, P-wave velocity, survey line azimuth angle, P-wave anisotropy parameter, P-S wave anisotropy parameter, fracture azimuth angle, and P-wave incident angle; A total travel time calculation model determination module, configured to obtain the total travel time calculation model by means of layer-by-layer inversion; the total travel time calculation model is a model of anisotropy parameters; the anisotropy parameters include P-wave anisotropy parameter, P-S wave anisotropy parameter, and S-wave anisotropy parameter; A first anisotropy parameter inversion module, configured to, based on the total travel time calculation model and the travel time information of the current formation at each survey line azimuth, use the least squares method to inversely solve the first objective function regarding the total travel time, and determine the optimal anisotropy parameters in the total travel time calculation model as the first inversion values of the anisotropy parameters; A first elastic parameter inversion module, configured to perform isotropic AVO inversion on the elastic parameters in the incident angle gather to obtain the first inversion values of the elastic parameters; the elastic parameters include P-wave velocity, S-wave velocity, and density; An anisotropic parameter secondary inversion module, configured to inversely solve a second objective function regarding the first inversion values of elastic parameters and the first inversion values of anisotropic parameters by using the Gauss-Newton method based on the first inversion values of anisotropic parameters and the first inversion values of elastic parameters, so as to obtain the secondary inversion values of anisotropic parameters; A prediction module, configured to predict fractures in the current formation according to the secondary inversion values of anisotropic parameters.

6. An electronic device, characterized in that, Comprising: One or more processors; A storage device storing one or more programs thereon; When the one or more programs are executed by the one or more processors, the one or more processors are caused to implement the method according to any one of claims 1 to 4.

7. A storage medium, characterized in that, A computer program is stored thereon, wherein when the computer program is executed by a processor, the method according to any one of claims 1 to 4 is implemented.

Images