Prestack crack prediction method and system under velocity anisotropy constraint

Abstract

The invention provides a pre-stack fracture prediction method and system under velocity anisotropy constraint, and relates to the technical field of seismic data processing, and the method comprises the steps: obtaining a seismic reflection coefficient approximation formula of HTI medium anisotropy parameters; according to the seismic reflection coefficient approximation formula, obtaining the relationship between the anisotropy parameter of the HTI medium and the velocity orientation, and through the inversion of the velocity orientation anisotropy, establishing an HIT medium seismic inversion objective function constrained by the velocity orientation anisotropy; and carrying out anisotropy inversion solution of amplitude difference on the HIT medium seismic inversion target function to obtain a target HTI medium anisotropy parameter representing anisotropy strength, thereby realizing amplitude difference prestack fracture prediction based on velocity and azimuth anisotropy constraint, and improving the fracture-vuggy carbonate rock prestack fracture prediction precision.

Description

Technical Field

[0001] This invention relates to the field of seismic data processing technology, and in particular to a method and system for predicting pre-stack cracks under velocity anisotropy constraints. Background Technology

[0002] Carbonate rocks are a common type of sedimentary rock. Compared to conventional sandstone and mudstone reservoirs, carbonate formation is influenced by tectonic and diagenetic processes, resulting in a complex and diverse range of carbonate reservoir types and significant differences in reservoir characteristics. Controlled by factors such as overlying strata pressure and strike-slip faults, carbonate reservoirs often develop numerous oriented vertical fractures along fault zones. The oriented arrangement of these fractures can be approximated as a transversely isotropic medium (HTI medium) with a horizontal axis of symmetry.

[0003] Seismic azimuth anisotropy inversion is crucial for the comprehensive interpretation of pre-stack wide-azimuth seismic data. The core of seismic inversion in complex carbonate fracture-vuggy reservoirs lies in how to stably extract anisotropic information from wide-azimuth seismic data. Existing anisotropic medium inversion methods lack consideration for the stability of simultaneous multi-parameter inversion in their inversion strategies and algorithms.

[0004] Furthermore, conventional pre-stack wide-azimuth fracture prediction techniques for carbonate reservoirs predict amplitude variations with azimuth and offset caused by fracture anisotropy. This method requires high aspect ratio and signal-to-noise ratio from pre-stack wide-azimuth seismic data. However, when both velocity and amplitude anisotropy are present in the wide-azimuth seismic data, the applicability of this method decreases, leading to a decline in the inversion accuracy of pre-stack fracture prediction. Summary of the Invention

[0005] The purpose of this invention is to provide a method and system for predicting pre-stack fractures under velocity anisotropy constraints, thereby improving the accuracy of pre-stack fracture prediction in fracture-cavity carbonate rocks.

[0006] To achieve the above objectives, the present invention provides the following technical solution:

[0007] In a first aspect, the present invention provides a method for predicting pre-stack cracks under velocity anisotropy constraints, comprising:

[0008] An approximate formula for obtaining the seismic reflection coefficient of the HTI medium anisotropy parameters;

[0009] The relationship between the anisotropy parameters and velocity azimuth of the HTI medium is obtained based on the approximate formula for the seismic reflection coefficient. Then, the objective function for seismic inversion of the HTI medium with velocity azimuth anisotropy constraint is established through the inversion of velocity azimuth anisotropy.

[0010] Anisotropic inversion of amplitude difference is performed on the seismic inversion objective function of the HIT medium to obtain the anisotropic parameters of the target HIT medium characterizing the anisotropic intensity, thereby realizing pre-stack crack prediction based on velocity azimuth anisotropy constraints.

[0011] Secondly, the present invention also provides a pre-stack crack prediction system under velocity anisotropy constraints, comprising:

[0012] Establish a unit to obtain an approximate formula for the seismic reflection coefficient of the HTI medium anisotropy parameters;

[0013] The inversion unit is used to obtain the relationship between the anisotropy parameters of the HTI medium and the velocity azimuth according to the approximate formula of the seismic reflection coefficient, and to establish the objective function of the HIT medium seismic inversion constrained by the velocity azimuth anisotropy through the inversion of velocity azimuth anisotropy.

[0014] The solution unit is used to perform amplitude difference anisotropy inversion solution on the seismic inversion objective function of the HIT medium, obtain the target HTI medium anisotropy parameters characterizing the anisotropy intensity, and realize pre-stack crack prediction based on velocity azimuth anisotropy constraints.

[0015] Thirdly, embodiments of the present invention also provide an electronic device, including a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program or instructions to implement the steps of the aforementioned method for predicting pre-stack cracks under velocity anisotropy constraints.

[0016] Fourthly, embodiments of the present invention also provide a computer storage medium storing a computer program or instructions, wherein when the computer program or instructions are executed by a processor, the steps of the aforementioned method for predicting pre-stack cracks under velocity anisotropy constraints are implemented.

[0017] Fifthly, embodiments of the present invention also provide a computer program product, including a computer program or instructions, which, when executed by a processor, implement the steps of the aforementioned pre-stack crack prediction method under velocity anisotropy constraints.

[0018] The technical effects and advantages of this invention are as follows: 1. This invention first derives the approximate formula for the reflection coefficient of HTI medium using the first-order scattering elastic wave steady-state method. Then, by establishing the relationship between anisotropic parameters and velocity, it establishes the inversion objective function constrained by velocity anisotropy. Then, based on the Bayesian framework, it performs anisotropic inversion of amplitude differences to solve the inversion objective functional. Finally, it realizes pre-stack fracture prediction based on velocity azimuth anisotropy constraints, further improving the fracture prediction accuracy of carbonate fracture-type reservoirs and providing technical support for efficient exploration and development of oil fields.

[0019] 2. This invention considers both velocity anisotropy and amplitude anisotropy caused by fracture anisotropy. The inversion method based on velocity azimuth anisotropy is relatively stable, but its resolution is low and it can only identify large reservoirs. On the other hand, the inversion method based on amplitude azimuth anisotropy has high resolution and is sensitive to the degree of medium anisotropy, but its noise resistance is relatively poor. Combining the advantages of both, this invention proposes a new fracture prediction technology. Under the constraint of velocity azimuth anisotropy, it combines the amplitude difference inversion method to carry out the inversion technology of HTI medium anisotropy parameters and fracture sensitive parameters, thereby improving the prediction accuracy, reliability and stability of pre-stack fracture prediction in fracture-vuggy carbonate reservoirs.

[0020] Other features and advantages of the invention will be set forth in the description which follows, and will be apparent in part from the description, or may be learned by practicing the invention. The objects and other advantages of the invention may be realized and obtained by means of the structures pointed out in the description, claims and drawings. Attached Figure Description

[0021] To more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0022] Figure 1 This invention provides a flowchart of a pre-stack crack prediction method under velocity anisotropy constraints. Figure 1 ;

[0023] Figure 2 This invention provides a flowchart of a pre-stack crack prediction method under velocity anisotropy constraints. Figure 2 ;

[0024] Figure 3 This is a schematic diagram of the structure of a pre-stack crack prediction system under velocity anisotropy constraint according to an embodiment of the present invention;

[0025] Figure 4 This is a schematic diagram of the structure of an electronic device according to an embodiment of the present invention;

[0026] Figure 5 This is a schematic diagram of the forward-modeling angle gathers of well A at different azimuth angles in an embodiment of the present invention;

[0027] Figure 6 This is a schematic diagram of the anisotropic parameter prediction results of well A under the velocity azimuth anisotropy constraint in an embodiment of the present invention.

[0028] Figure 7 This is a schematic diagram of seismic profiles passing through well A at different azimuth angles in an embodiment of the present invention;

[0029] Figure 8 This is a schematic diagram of the fracture anisotropy parameter inversion profile through well A under velocity azimuth anisotropy constraint in an embodiment of the present invention.

[0030] Figure 9 This is a schematic diagram of a plane slice for inverting crack anisotropic parameters under velocity-azimuth anisotropy constraints in an embodiment of the present invention. Detailed Implementation

[0031] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0032] Due to the characteristics of hydrochloric acid reservoirs, such as large burial depth, dense rocks, complex pore structure, strong heterogeneity and anisotropy, the relationship between reservoir anisotropy parameters and seismic response characteristics is complex.

[0033] To clarify the seismic response characteristics of fractured carbonate reservoirs (HTI media) with directional arrangement, this embodiment of the invention obtains the correspondence between the anisotropic parameters of fractured carbonate reservoirs and seismic response characteristics by acquiring the approximate formula of the reflection coefficient of the anisotropic parameters of HTI media, that is, to obtain the actual observation data. This process is the solution process of the forward problem.

[0034] Therefore, this invention uses observational data from a wide azimuth to inversely deduce the anisotropy parameters of fractured-vuggy carbonate reservoirs with HTI characteristics, a process that is an inverse problem solution.

[0035] For fractured-vuggy carbonate reservoirs with strong heterogeneity and HTI characteristics, conventional post-stack fracture prediction methods are not accurate enough to meet the needs of fine fracture prediction.

[0036] To improve the accuracy of pre-stack fracture prediction in fracture-vuggy carbonate rocks and reduce the ambiguity of inversion results, this invention discloses a pre-stack fracture prediction method under velocity anisotropy constraints, such as... Figure 1 and Figure 2 As shown, it includes the following steps:

[0037] Step S1: For fractured-vuggy carbonate HTI media, obtain the approximate formula for the seismic reflection coefficient of the anisotropic parameters of HTI media (i.e., the fracture parameters of HTI media) to establish a quantitative relationship between the seismic reflection coefficient and the anisotropic parameters of HTI media.

[0038] Step S2: Use the approximate formula for the seismic reflection coefficient to obtain the relationship between the anisotropic parameters of the HTI medium and the velocity azimuth, and establish the objective function for the seismic inversion of the HTI medium with velocity azimuth anisotropy constraint through the inversion of velocity azimuth anisotropy, thereby realizing the velocity anisotropy relationship.

[0039] Step S3: Perform amplitude difference anisotropy inversion solution on the seismic inversion objective function of the HIT medium to obtain the target anisotropy parameters characterizing the anisotropy intensity, and realize pre-stack crack prediction based on velocity azimuth anisotropy constraint.

[0040] In some specific embodiments, step S1: For fractured-vuggy HTI carbonate media, obtain an approximate formula for the seismic reflection coefficient of the HTI media anisotropy parameters, so as to establish a quantitative relationship between the seismic reflection coefficient and the HTI media anisotropy parameters; the specific operation is as follows:

[0041] In order to solve the inverse problem, the forward problem is described first. Since constructing the forward operator is the key to the anisotropic inversion of HTI media, this invention first obtains the approximate formula of the reflection coefficient characterized by the anisotropic parameters, so as to facilitate the subsequent construction of the forward operator.

[0042] Therefore, this invention focuses on fractured-vuggy carbonate rock HTI media containing a single set of vertical fractures. Based on Schoenberg linear slip theory, the relationship between the scattering function and the reflection characteristic equation of the HTI medium is obtained using the first-order scattering elastic wave steady-state method:

[0043]

[0044] In the formula, R pp (θ) represents the approximate formula for the reflection coefficient of PP waves, ρ b S represents the background medium density. PP (r0) represents the scattering function value at point r = r0, where r represents the position of the scattering function and r0 represents the initial position of the scattering function.

[0045] Wherein, the scattering function S PP (r0) is:

[0046]

[0047] In the formula, Δρ represents the density perturbation term, and ΔC mnξ represents the perturbation term of the elastic matrix of the HTI medium. pp This represents the contribution of the incident PP wave to the scattering function. This represents the contribution of the scattered PP wave to the scattering function; the relationship between the subscripts m and n and i, j, k, and l can be expressed as:

[0048] In the formula, δ ij and δ kl Both refer to the Kronecker function;

[0049] ξ pp and It can be represented by the slowness p and g of the incident PP wave and the scattered PP wave, and the polarization vectors p′ and g′, that is:

[0050]

[0051] In the formula, This represents the polarization vector of the incident PP wave along the x-direction. This indicates the slowness of the incident PP wave along the y-direction. This represents the polarization vector of the incident PP wave along the y-direction. Let r represent the slowness of the incident PP wave along the z-direction, r represent the position of the scattering function, and r0 represent the initial position of the scattering function.

[0052]

[0053] Where φ represents porosity, Represents the imaginary part along the x-direction. This represents the imaginary part along the y-direction. α represents the imaginary part of the table along the z-direction. b Represents the longitudinal wave velocity of the background medium;

[0054] therefore, Among them, t i This represents the sampling time of the i-th sampling point in the observed seismic record, where i represents the sequence number (here, the sequence number of the sampling point); t′ i Let represent the propagation time of the i-th sampling point, r represent the position of the scattering function, and r0 represent the initial position of the scattering function; and

[0055] η 11 =(sin 4 θcos 4 φ) / vp 2 η 12 =(sin 4 θsin 2 φcos 2 φ) / vp 2 ,

[0056] or 13 =(sin 2 θcos 2 θcos 2 f) / vp 2 ,or 14 =2(sin 3 θcosθsinφcos 2 f) / vp 2 ,

[0057] or 15 =-2(sin 3 θcosθcos 3 f) / vp 2 ,or 16 =2(sin 4 θsinφcos 3 f) / vp 2 ,

[0058] or 21 =(sin 4 θsin 2 φcos 2 f) / vp 2 ,or 22 =(sin 4 θsin 4 f) / vp 2 ,

[0059] or 23 =(sin 2 θcos 2 θsin 2 f) / vp 2 ,or 24 =-2(sin 3 θcosθsin 3 f) / vp 2 ,

[0060] or 25 =-2(sin 3 θcosθsin 2 φcosφ) / vp 2 ,or 26 =2(sin 4 θsin 3 φcosφ) / vp 2 ,

[0061] or 31 =(sin 2 θcos 2 θcos 2 f) / vp 2 ,or32 =(sin 2 θcos 2 θsin 2 f) / vp 2 ,

[0062] or 33 =(cos 4 i) / vp 2 ,or 34 =2(sinθcos 3 θsinφ) / vp 2 ,

[0063] or 35 =-2(sinθcos 3 θcosφ) / vp 2 ,or 36 =2(sin 2 θcos 2 θsinφcosφ) / vp 2 ,

[0064] or 41 =-2(sin 3 θcosθsinφcos 2 f) / vp 2 ,or 42 =2(sin 3 θcosθsin 3 f) / vp 2 ,

[0065] or 43 =-2(sinθcos 3 θsinφ) / vp 2 ,or 44 =-4(sin 2 θcos 2 θsin 2 f) / vp 2 ,

[0066] or 45 =-4(sin 2 θcos 2 θsinφcosφ) / vp 2 ,or 46 =-4(sin 3 θcosθsin 2 φcosφ) / vp 2 ,

[0067] or 51 =2(sin 3 θcosθcos 3 f) / vp2 ,or 52 =2(sin 3 θcosθsin 2 φcosφ) / vp 2 ,

[0068] or 53 =2(sinθcos 3 θcosφ) / vp 2 ,or 54 =-4(sin 2 θcos 2 θsinφcosφ) / vp 2 ,

[0069] or 55 =-4(sin 2 θcos 2 θcos 2 f) / vp 2 ,or 56 =-4(sin 3 θcosθsinφcos 2 f) / vp 2 ,

[0070] or 61 =2(sin 4 θsinφcos 3 f) / vp 2 ,or 62 =2(sin 4 θsin 3 φcosφ) / vp 2 ,

[0071] or 63 =2(sin 2 θcos 2 θsinφcosφ) / vp 2 ,or 64 =4(sin 3 θcosθsin 2 φcosφ) / vp 2 ,

[0072] or 65 =4(sin 3 θcosθsinφcos 2 f) / vp 2 ,or 66 =4(sin 4 θsin 2 φcos 2 f) / vp 2 。

[0073] Where φ represents porosity and vp represents the longitudinal wave velocity of the background medium.

[0074] Because of the scattering function S in formula (2) PP Since there is a perturbation term in (r0), the stiffness matrix of the HTI medium can be expressed in perturbation form based on the perturbation term in the scattering function. That is, the perturbation stiffness matrix of the HTI medium is obtained as follows:

[0075]

[0076] in,

[0077]

[0078] ΔC 44 =-qΔφμ m +Δμ m

[0079] ΔC 55 =-qΔφμ m -μ m Δδ T +Δμ m ,

[0080] In the formula, ΔC HTI The symbol Δ represents the disturbance stiffness of the HTI medium, φ represents porosity, Δφ represents the porosity disturbance of the rock, P represents the effective pressure of the rock, ΔP represents the effective pressure disturbance of the rock, and K represents the perturbation stiffness of the rock. m ΔK represents the bulk modulus of the rock matrix. m δ represents the disturbance of the bulk modulus of the rock matrix. N With δ T Both represent the rock physical parameters of the fracture in the Schoenberg linear slip model, namely the fracture normal weakness and fracture tangential weakness; Δδ N With Δδ T μ represents the weak perturbation in the normal direction and the weak perturbation in the tangential direction of the crack, respectively. m Δμ represents the shear modulus of the rock matrix. m This represents the amount of shear modulus disturbance in the rock matrix; x represents the ratio of the Lamé parameter to the bulk modulus of the rock. M = λ + 2μ, where M represents the bulk modulus of the rock, λ represents the Lamé parameter, and μ represents the shear modulus of the rock.

[0081] Substituting the perturbation stiffness matrix of the HTI medium into the relationship between the scattering function and the reflection characteristic equation of the HTI medium, an approximate formula for the reflection coefficient of the crack parameters in the HTI medium is obtained, which is derived from the bulk modulus K of the rock matrix. m Shear modulus μ of the rock matrix mThe rock's density ρ, porosity φ, effective pressure P, and fracture normal weakness δ N and crack tangential weakness δ T Approximate formula for the seismic reflection coefficient for:

[0082]

[0083] in,

[0084]

[0085] In the formula, The reflection coefficient of a fractured reservoir varies with the incident angle θ and azimuth angle. The change in θ represents the angle of incidence, i.e., the angle of incidence of the seismic wave; The azimuth angle represents the angle between the azimuth angle of the observation line and the axis of symmetry of the crack. φ represents the observable azimuth angle, φ sym Indicates the symmetrical azimuth angle of the crack; The bulk modulus K of the rock matrix m The coefficient K that varies with the incident angle θ m0 This represents the mean bulk modulus of the rock matrix. The shear modulus μ of the rock matrix m The coefficient a that varies with the incident angle θ ρ (θ) represents the coefficient by which the density ρ of the rock varies with the incident angle θ, Δρ represents the disturbance of the rock density, ρ0 represents the mean density of the rock, and a φ (θ) represents the coefficient by which the porosity φ of the rock varies with the incident angle θ, φ0 represents the mean porosity of the rock, and a p (θ) represents the coefficient of variation of the effective pressure P of the rock with the incident angle θ, and P0 represents the mean effective pressure of the rock. This indicates that the crack normal weakness varies with the incident angle θ and azimuth angle. The coefficient of change This indicates that the tangential weakness of the crack varies with the incident angle θ and the azimuth angle. The coefficient of variation, g1, represents the ratio of the bulk modulus of the rock matrix to the longitudinal wave modulus of the rock, i.e., g1 = K. m / M,K m g1 represents the bulk modulus of the rock matrix; g2 represents the ratio of the shear modulus of the rock matrix to the longitudinal wave modulus of the rock, i.e., g2 = μ. m / M,μ m g3 represents the shear modulus of the rock matrix; g3 represents the ratio of the rock shear modulus to the longitudinal wave modulus. Where M and μ represent the longitudinal wave modulus and shear modulus of the rock, respectively, and g1, g2 and g3 can be estimated using well logging data.

[0086] Therefore, formula (7) gives an approximate formula for the seismic reflection coefficient that directly characterizes the anisotropy parameters of the HTI medium, thus establishing the solution process for the forward problem.

[0087] In this embodiment of the invention, step S1 is used to establish the relationship between the seismic reflection coefficient and the anisotropy parameters of the HTI medium, thereby achieving amplitude anisotropy.

[0088] Conventional inversion methods based on amplitude azimuth anisotropy offer advantages such as high resolution and sensitivity to the degree of medium anisotropy, but their resistance to noise interference is relatively poor, requiring high-quality seismic data. Inversion methods based on velocity azimuth anisotropy are more stable, but their resolution is lower, and they can only identify large reservoirs.

[0089] Current conventional inversion methods based on amplitude azimuth anisotropy lack consideration for velocity azimuth anisotropy, necessitating the development of more effective stable inversion methods for anisotropic media parameters. The propagation velocity of seismic waves in azimuthally anisotropic media changes with the propagation direction. Using the optimal imaging reference velocity, after normal time difference correction, the reflection phase axis on the gather should tend to flatten. Conversely, the residual time difference is caused by azimuth anisotropy. In short, azimuthally anisotropic residual time difference refers to the perturbation time difference after the reflection event is flattened, based on which velocity azimuth anisotropy information can be obtained.

[0090] Therefore, in some specific embodiments, step S2 involves obtaining the relationship between anisotropy parameters and velocity azimuth using the approximate formula for the seismic reflection coefficient, and establishing a seismic inversion objective function for HIT media constrained by velocity azimuth anisotropy through the inversion of velocity azimuth anisotropy; the specific operation is as follows:

[0091] Step S21: Based on the aforementioned approximate formula for the reflection coefficient, obtain the fracture rock physical parameters δ of the Schoenberg linear slip model. N and δ T ,include:

[0092] Based on the aforementioned approximate formula for the reflection coefficient, the seismic P-wave phase velocity of the HIT medium, characterized by the Thomsen anisotropy parameter, is as follows:

[0093]

[0094] In the formula, Let θ represent the incident angle of the seismic wave and azimuth angle be θ. Seismic P-wave phase velocity in HIT medium The value represents the phase velocity of the seismic P-wave, and θ represents the angle of incidence, i.e., the angle of incidence of the seismic wave. The azimuth angle represents the angle between the azimuth angle of the observation line and the axis of symmetry of the crack. φ Observe Indicates the observable azimuth angle, φ sym Indicates the symmetrical azimuth angle of the crack; V is o represents the isotropic velocity of the seismic P-wave, ε (V) δ (V) and γ (V) All represent the Thomsen anisotropy parameters in the VTI background.

[0095] By utilizing the relationship between the Schoenberg linear sliding model and the Hudson thin coin-shaped crack model, the Thomsen anisotropy parameter ε is obtained through inversion. (V) δ (V) and γ (V) Rock physical parameters δ of fractures compared with the Schoenberg linear slip model N and δ T The relationship between them is:

[0096]

[0097] In the formula, ε (V) δ (V) and γ (V) Both represent the Thomsen anisotropy parameter, and g represents the square of the ratio of seismic S-wave to P-wave velocity under isotropic background. V P0 V represents the phase velocity of the seismic P-wave in the HIT medium. S δ represents the transverse wave velocity of the rock. N and δ T Both represent the rock physics parameters of the Schoenberg fracture model, and are also two parameters included in the HTI medium anisotropy parameters.

[0098] Based on Thomsen anisotropy parameter ε (V) δ (V) and γ (V) Rock physical parameters δ of fractures compared with the Schoenberg linear slip model N and δ T Based on the relationship between these parameters and the seismic P-wave phase velocity characterizing the HIT medium, the fracture rock physical parameters of the Schoenberg linear slip model are obtained as follows:

[0099]

[0100] In the formula, δ N and δ T Both represent the rock physical parameters of the fracture in the Schoenberg linear slip model, namely the weakness in the normal direction and the weakness in the tangential direction of the fracture; ε(V) δ (V) and γ (V) Both represent the Thomsen anisotropy parameter, and g represents the square of the ratio of seismic S-wave to P-wave velocity under isotropic background.

[0101] Step S22: Using the anisotropic parameters of the HTI medium to be inverted as the parameters to be inverted, the velocity orientation of the HTI medium anisotropic parameters is inverted to obtain the HTI medium anisotropic parameter matrix.

[0102] Based on the Schoenberg linear slip model, the rock physical parameters δ of the fracture are... N and δ T And the anisotropic parameter matrix of the HTI medium, to obtain the velocity orientation and anisotropic constraints:

[0103] Π vel =λ i (η i -l i m i ) Τ (η i -l i m i (11),

[0104] In the formula, λ i δ N and δ T constraint coefficient, in, δ N constraint coefficient, δ T constraint coefficient; η i η represents the scaling factor between the parameters to be inverted and their initial values. i =1 / 2×ln(m) i / m i0 ), m i0 The initial value of the anisotropy parameter matrix obtained from velocity azimuth inversion is represented by l. i This represents the integral operator in a mathematical formula. t i This represents the sampling time of the i-th sampling point in the observed seismic record, where i represents the sequence number (here, the sequence number of the sampling point); t0 represents the initial sampling time; m i This represents the anisotropic parameter matrix obtained from velocity azimuth inversion.

[0105] Step S23: Based on the velocity azimuth and anisotropy constraints, establish the seismic inversion objective function for the HIT medium as follows:

[0106]

[0107] In the formula, Let ΔS represent the seismic inversion objective function for the HIT medium, ΔG′ represent the perturbation of the seismic record matrix, and ΔG′ represent the forward modeling operator matrix for azimuth amplitude difference. This represents the reflection coefficient matrix of the parameters to be inverted after decorrelation. The variance of the random noise is represented by K, the total number of seismic traces involved in the inversion is represented by i, and the trace number is represented here. i This represents the anisotropic parameter matrix obtained from velocity azimuth inversion. Π represents the variance of the parameter to be inverted. vel This indicates the velocity orientation and anisotropy constraints.

[0108] In some specific embodiments, step S3: The amplitude difference anisotropy is inverted and solved for the seismic inversion objective function of the HIT medium to obtain the target HTI medium anisotropy parameters characterizing the anisotropy intensity, thereby realizing pre-stack crack prediction based on velocity azimuth anisotropy constraints. The specific operation is as follows:

[0109] In order to optimize the objective function of the seismic inversion of the HIT medium, this embodiment of the invention combines a Bayesian inversion framework and uses the iterative reweighted least squares algorithm to solve it.

[0110] That is, step S31: Under the Bayesian inversion framework, the anisotropic parameters to be inverted are used as the parameters to be inverted. Assuming that the parameters to be inverted and the noise in the observed seismic records follow Cauchy and Gaussian distributions respectively, the posterior probability density distribution function based on the parameters to be inverted and the observed seismic records is obtained, specifically including:

[0111] Bayes' theorem is a theorem concerning the conditional or marginal probabilities of random events A and B. Generally, the probability of event A given that event B has occurred is different from the probability of event B given that event A has occurred; however, there is a definite relationship between the two, and Bayes' theorem states this relationship. One application of Bayes' formula is to deduce a fourth probability function from three known probability functions.

[0112]

[0113] In the formula, P(A|B) represents the conditional probability of A after B has occurred, which is called the posterior probability of A; i represents the index, j represents the index, P(A) represents the prior probability or marginal probability of A, which does not consider any factors related to B; similarly, P(B|A) represents the conditional probability of B after A has occurred, which is called the posterior probability of A; P(B) represents the prior probability or marginal probability of B, which is also a standardized constant.

[0114] Since the Cauchy distribution is a long-tailed distribution function, it has high resolution in an algorithmic sense, can protect weak reflections such as anisotropic parameters, and increase the stability of anisotropic parameter inversion.

[0115] Therefore, in this embodiment of the invention, it is assumed that the parameters to be inverted follow a Cauchy distribution, and the probability distribution of the parameters to be inverted is as follows:

[0116]

[0117] In the formula, p Cauchy (m) represents the probability distribution of the parameters to be inverted, M represents the number of sampling points in the observed seismic records, and i represents the sequence number (here, the sequence number of the sampling point). R represents the variance of the parameter vector m to be inverted; i Let R represent the reflection coefficient corresponding to the i-th sampling point, and let R represent the reflection coefficient of the parameter vector m to be inverted.

[0118] Assuming the noise in the observed seismic records follows a Gaussian distribution, the likelihood function p(S|m) of the observed seismic records and the parameter matrix to be inverted is determined as follows:

[0119]

[0120] In the formula, p(S|m) represents the likelihood function of the observed seismic record and the parameter matrix to be inverted, and σ k This represents the noise distribution in the observed seismic record, where S represents the observed seismic record containing noise, i.e., the observed seismic record. denoted by , where represents the variance of the noise distribution in the observed seismic record, G represents the joint matrix of the wavelet matrix and the coefficient matrix of the parameter to be inverted, and m represents the parameter matrix to be inverted;

[0121] Based on the probability distribution of the parameters to be inverted and the likelihood function, the posterior probability density distribution function based on the parameters to be inverted and the noisy seismic record is obtained as follows:

[0122]

[0123] In the formula, p(S|m) represents the likelihood function of the observed seismic record and the parameter matrix to be inverted, W represents the number of sampling points in the observed seismic record, and i represents the index (here, the index of the sampling point); R i This represents the reflection coefficient corresponding to the i-th sampling point. The variance of the parameter vector m to be inverted is represented by S, and the observed seismic record containing noise is represented by the observed seismic record. Let represent the variance of the noise distribution in the observed seismic record, G represent the joint matrix of the wavelet matrix and the coefficient matrix of the parameter to be inverted, and m represent the parameter matrix to be inverted.

[0124] Step S32: Using the iterative reweighted least squares algorithm and the posterior probability density distribution function, solve the objective function of the seismic inversion of the HIT medium to obtain the parameter matrix m to be inverted.

[0125] Step S33: Based on the simplified approximate formula for the reflection coefficient, obtain the forward modeling equation of the HTI medium observation seismic record characterized by the anisotropy parameters of the HTI medium, and obtain the observation seismic record matrix S according to the forward modeling equation of the HTI medium observation seismic record; including the following steps:

[0126] Step S331: Simplify the approximate formula for the reflection coefficient, and obtain the forward modeling equation of the HTI medium seismic record characterized by the anisotropy parameters of the HTI medium based on the simplified approximate formula for the reflection coefficient. That is, obtain the forward modeling equation of the HTI medium seismic record characterized by the parameters to be inverted, specifically including:

[0127] The relative differences of the background medium parameters to be inverted in the approximate formula for the reflection coefficient are expressed in the form of the difference of their natural logarithms, i.e. Given Δρ / ρ0≈Δ(lnρ), Δφφ0≈Δ(lnφ), and ΔP / P0≈Δ(lnP), by removing the fractional terms from the aforementioned approximate formula for the reflection coefficient, we obtain the simplified approximate formula for the reflection coefficient as follows:

[0128]

[0129] In the formula, R pp (θ,φ) represents the reflection coefficient of a fractured reservoir as a function of the incident angle θ and azimuth angle. Changes, The bulk modulus K of the rock matrix m The coefficient that varies with the incident angle θ. a is a coefficient representing the variation of the shear modulus of the rock matrix with the incident angle θ. ρ (θ) represents the coefficient that the density ρ of the rock varies with the incident angle θ, a φ (θ) represents the coefficient that the porosity φ of the rock varies with the incident angle θ, a p (θ) represents the coefficient that represents the effective pressure P of the rock as a function of the incident angle θ. Indicates the crack normal weakness δ N With incident angle θ and azimuth angle The coefficient of change Indicates the tangential weakness δ of the crack T With incident angle θ and azimuth angle The coefficient of change The bulk modulus K of the rock matrix m Contribution to the reflection coefficient K m Indicates the bulk modulus of the rock matrix; The shear modulus μ of the rock matrix m Contribution to the reflection coefficient μ m R represents the shear modulus of the rock matrix. ρ R represents the contribution of rock density ρ to the reflectance coefficient. ρ =Δln(ρ), where ρ represents the density of the rock; R φ R represents the contribution of rock porosity φ to the reflection coefficient. φ =Δln(φ), where φ represents the porosity of the rock; R p R represents the contribution of the effective pressure P of the rock to the reflection coefficient. P =Δln(P), where P represents the effective pressure; Indicates the crack normal weakness δ N The contribution to the reflection coefficient, Indicates the tangential weakness δ of the crack T The contribution to the reflection coefficient, δ N With δ T Both represent the rock physical parameters of the fracture in the Schoenberg linear slip model, namely the fracture normal weakness and fracture tangential weakness.

[0130] By convolving the simplified reflection coefficient approximation formula with the seismic wavelet, the forward modeling equation for HTI medium seismic records, characterized by the anisotropic parameters of the HTI medium, is obtained as follows:

[0131]

[0132] In the formula, S represents the observed seismic record containing noise. This represents seismic wavelets from different directions and angles.

[0133] Step S331: Utilize isotropic background seismic records S is o and anisotropic seismic records The observed seismic records in the forward modeling equations of the HTI medium seismic records are represented in the form of a sum, and then converted into matrix form to obtain the observed seismic record matrix S, specifically including:

[0134] Using isotropic background seismic records S iso With anisotropic earthquake records The observed seismic records in the forward modeling equations of the HTI medium seismic records can be represented in the form of a sum as follows:

[0135]

[0136] In the formula, S represents the observed seismic record containing noise, S iso This represents isotropic background earthquake records. Represents anisotropic earthquake records;

[0137] The observed seismic record matrix is ​​as follows:

[0138]

[0139] In the formula, [S] represents the observed seismic record matrix, G represents the joint matrix of the wavelet matrix and the coefficient matrix of the parameters to be inverted, and G = Wa(θ); G iso The joint coefficient matrix of the isotropic background parameters is represented. Let m represent the joint coefficient matrix of the parameters to be inverted. iso This represents the parameter matrix to be inverted against an isotropic background. The bulk modulus K of the rock matrix m Contribution to the reflection coefficient The shear modulus μ of the rock matrix m The contribution of R to the reflection coefficient ρ R represents the contribution of rock density ρ to the reflectance coefficient. φ R represents the contribution of rock porosity φ to the reflection coefficient. p This represents the contribution of the effective pressure P of the rock to the reflection coefficient. This represents the anisotropic parameter matrix to be inverted. Indicates the crack normal weakness δ N The contribution to the reflection coefficient, Indicates the tangential weakness δ of the crack T The contribution to the reflection coefficient.

[0140] Step S34: Based on Bayesian theory, the observed seismic record matrix S is transformed using the anisotropic medium AVAZ inversion method, the decorrelated coefficient matrix G′, and the decorrelated parameter matrix m′ to be inverted, to obtain the decorrelated observed seismic record matrix S′. Specifically, this includes:

[0141] Based on Bayesian theory, the isotropic forward modeling operator G is obtained using the AVAZ inversion method for anisotropic media. iso and anisotropic orthogonal operators

[0142] When the azimuth seismic gather has N incident angles, M azimuth angles, and F reflecting interfaces, the isotropic forward modeling operator G... iso for:

[0143]

[0144] In the formula, [G iso ] NMF×5F G represents the isotropic forward modeling operator when the azimuth seismic gather has N incident angles, M azimuth angles, and F reflecting interfaces. iso W represents the wavelet matrix. This represents the matrix of bulk modulus coefficients of the rock matrix at N incident angles. Let A represent the matrix of rock matrix shear modulus coefficients for N incident angles. ρ (θ) represents the density coefficient matrix of the rock at N incident angles, A φ (θ) represents the porosity coefficient matrix of the rock at N incident angles, A P (θ) represents the effective pressure coefficient matrix of the rock at N incident angles;

[0145] Anisotropic forward modeling operator when the azimuth seismic gather has N incident angles, M azimuth angles, and F reflecting interfaces. for:

[0146]

[0147] In the formula, This represents the anisotropic forward modeling operator when the azimuth seismic gather has N incident angles, M azimuth angles, and F reflecting interfaces. W represents the wavelet matrix. This represents the matrix of crack normal weakness coefficients for N incident angles and M azimuth angles. This represents the crack tangential weakness coefficient matrix for N incident angles and M azimuth angles.

[0148] To improve the stability of inversion predictions, the covariance matrix C is used. c For the isotropic forward operator G iso Anisotropic orthogonal operators The parameter matrix m to be inverted is then subjected to decorrelation processing to obtain the decorrelation coefficient matrix G′ and the parameter matrix m′ to be inverted as follows:

[0149]

[0150] In the formula, G i ′ so G represents the de-correlated isotropic forward modeling operator. iso This represents an isotropic forward modeling operator. This indicates the de-correlated anisotropic forward modeling operator. Denotes the anisotropic forward modeling operator, m i ′ so m represents the parameter matrix to be inverted from the isotropic background to be removed. isoThis represents the parameter matrix to be inverted against an isotropic background. This represents the anisotropic parameter matrix to be inverted to remove correlation. represents the anisotropic parameter matrix to be inverted, and u represents the eigenvector;

[0151] The seismic record matrix is ​​transformed using the decorrelation coefficient matrix G′ and the decorrelation parameter matrix m′ to obtain the decorrelation observed seismic record matrix S′:

[0152]

[0153] In the formula, [S′] represents the decorrelated observation seismic record matrix, and G i ′ so This indicates the de-correlated isotropic forward modeling operator. Denotes the decorrelated anisotropic forward modeling operator, m i ′ so This represents the parameter matrix to be inverted against the isotropic background of the HTI medium. This represents the inversion parameter matrix for the anisotropic HTI medium, which is related to the dereferenced parameters.

[0154] Step S35: Based on the decorrelation coefficient matrix G′ and the decorrelation observed seismic record matrix S′, the objective function for seismic inversion of the HIT medium is optimized to obtain the objective function for minimizing azimuth amplitude difference:

[0155]

[0156] In the formula λ represents the minimum azimuth amplitude difference coefficient, ΔG' represents the azimuth amplitude difference forward modeling operator matrix, and λ represents the minimum azimuth amplitude difference coefficient. c Q represents the ratio of the variance of the noise distribution in the seismic observation record to the variance of the parameters to be inverted. c This represents the diagonal matrix of variance coefficients of the inversion parameters. δ N constraint coefficient, δ N The probability distribution function, δ T constraint coefficient, δ T The probability distribution function, ΔS represents the contribution of the parameters to be inverted to the observed seismic record, and ΔS represents the composite seismic record of azimuth amplitude difference, that is, the amount of disturbance to the observed seismic record.

[0157] in, σ represents the variance of the noise distribution in earthquake observation records. mThis represents the variance of the parameters to be inverted, i.e., the variance of the anisotropic parameters to be inverted; m n σ represents the weighting coefficient of the parameters to be inverted. R' This represents the variance of the parameter vector m to be inverted;

[0158] Furthermore, the iterative reweighted least squares algorithm is used to further solve the objective function of small-azimuth amplitude difference to obtain the HTI medium anisotropy parameter δ, which characterizes the anisotropy intensity. N and δ T .

[0159] Based on the same inventive concept, embodiments of the present invention also provide a pre-stack crack prediction system under velocity anisotropy constraints, such as... Figure 3 As shown, it includes:

[0160] Establish a unit to obtain an approximate formula for the seismic reflection coefficient of the HTI medium anisotropy parameters;

[0161] The inversion unit is used to obtain the relationship between the anisotropy parameters of the HTI medium and the velocity azimuth according to the approximate formula of the seismic reflection coefficient, and to establish the objective function of the HIT medium seismic inversion constrained by the velocity azimuth anisotropy through the inversion of velocity azimuth anisotropy.

[0162] The solution unit is used to perform amplitude difference anisotropy inversion solution on the seismic inversion objective function of the HIT medium, obtain the target HTI medium anisotropy parameters characterizing the anisotropy intensity, and realize pre-stack crack prediction based on velocity azimuth anisotropy constraints.

[0163] Regarding the system in the above embodiments, the specific manner in which each unit module performs operations has been described in detail in the embodiments related to the method, and will not be elaborated here.

[0164] Based on the same inventive concept, embodiments of the present invention also provide an electronic device, the structure of which is as follows: Figure 4 As shown, it includes a memory, a processor, and a computer program stored in the memory. The processor executes the computer program or instructions to implement the steps of the aforementioned pre-stack crack prediction method under velocity anisotropy constraints.

[0165] Based on the same inventive concept, embodiments of the present invention also provide a computer storage medium storing a computer program or instructions, which, when executed by a processor, implement the steps of the aforementioned pre-stack crack prediction method under velocity anisotropy constraints.

[0166] Based on the same inventive concept, embodiments of the present invention also provide a computer program product, including a computer program or instructions. When the computer program or instructions are executed by a processor, they implement the steps of the aforementioned pre-stack crack prediction method under velocity anisotropy constraints. The computer program or instructions are stored in a computer storage medium. When the processor of an electronic device reads the computer program or instructions from the computer storage medium and executes the computer program or instructions, the electronic device performs the method described above in the embodiments of the present invention.

[0167] To ensure the rationality of the method in this invention, the approximate reflection coefficient formula (7) derived in step 1 is first used to perform a convolution model with the Ricker wavelet to obtain the synthetic seismic record angle gather as shown in the appendix. Figure 5 As shown; Appendix Figure 5 Angle gather records for well A at five azimuth angles of 30°, 60°, 90°, 120° and 150° are given respectively. As can be seen from the figure, there are some differences in gather amplitude at different azimuth angles. This invention uses the amplitude differences between different azimuth angles to carry out pre-stack fracture prediction.

[0168] To verify the reliability of the method in this invention, pre-stack crack prediction was performed on the forward-modeled gathers at the five different azimuth angles described above using the method provided in the embodiments of this invention. The prediction results are attached. Figure 6 As shown, from Figure 6 As can be seen from the figure, the dashed line represents the initial model of anisotropic parameters obtained by velocity anisotropy constraint (i.e., the prior information in Formula 13; this application requires adding an initial model as a constraint to solve the inverse problem), and the solid line represents the comparison of anisotropic parameters between the inversion result and the measured results on the well. The comparison shows that the anisotropic parameters predicted by the method of this invention are in good agreement with the measured results on the well, thus confirming the effectiveness of the technical method from the two-dimensional model.

[0169] To further verify the reliability of the method in this invention, the pre-stack fracture prediction method was applied to an actual work area of ​​a three-dimensional carbonate fracture-vuggy reservoir. Figure 7 The seismic data from well A at azimuth angles of 15°, 45°, 75°, 105°, 135°, and 165° are compared. Figure 7 As can be seen, seismic data at different azimuth angles exhibit amplitude differences across the profile. Using seismic data from different azimuth angles, pre-stack fracture prediction with velocity anisotropy constraints was performed. The predicted fracture anisotropy parameters through well A are shown in the attached figure. Figure 8 As shown, from Figure 8As can be seen from the data, the fracture anisotropy parameter inversion results have good continuity and high resolution. The high-value anomaly responses are well-matched with the fracture development regions in the well logging interpretation. Therefore, anisotropy parameters can be used to indicate fracture development regions. A planar plot of the fracture anisotropy parameter inversion results is attached. Figure 9 As shown in the sectional slice results, the fracture anisotropy parameters at the well location in the work area are all high-value anomalies, which are in good agreement with the well logging interpretation results, further confirming the rationality of the method provided by the present invention in the actual work area.

[0170] To improve the reliability of multi-parameter direct inversion, this invention, based on Bayesian theory, employs an azimuth amplitude difference inversion strategy to conduct azimuth amplitude difference inversion of fractured-vuggy carbonate reservoirs in HTI media. Driven by pre-stack wide-azimuth seismic data, it studies the inversion methods for seismic anisotropy parameters and fracture-sensitive parameters. Since current inversion methods based on amplitude azimuth anisotropy lack consideration for velocity azimuth anisotropy, there is an urgent need to develop more effective stable inversion methods for anisotropic media parameters. Therefore, to further improve the stability of anisotropic parameter inversion, this invention, based on velocity azimuth anisotropy characteristics and combined with velocity azimuth anisotropy constraints, proposes a new amplitude difference inversion method to conduct inversion of HTI media anisotropic parameters and fracture-sensitive parameters. Tests using two-dimensional models and actual data demonstrate that the inversion results of the proposed method have a high degree of agreement with well logging interpretation results, proving that this inversion method has good application prospects in anisotropic media inversion.

[0171] Finally, it should be noted that the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for predicting pre-stack cracks under velocity anisotropy constraints, characterized in that, include: An approximate formula for obtaining the seismic reflection coefficient of the HTI medium anisotropy parameters; The relationship between the anisotropy parameters and velocity azimuth of the HTI medium is obtained based on the approximate formula for the seismic reflection coefficient. Then, the objective function for seismic inversion of the HTI medium with velocity azimuth anisotropy constraint is established through the inversion of velocity azimuth anisotropy. Anisotropic inversion of amplitude difference is performed on the seismic inversion objective function of the HIT medium to obtain the anisotropic parameters of the target HIT medium characterizing the anisotropic intensity, thereby realizing pre-stack crack prediction based on velocity azimuth anisotropy constraints.

2. The method for predicting pre-stack cracks under velocity anisotropy constraints according to claim 1, characterized in that, The approximate formula for obtaining the reflection coefficient of the HTI medium anisotropy parameters includes: The relationship between the scattering function and the HTI medium reflection characteristic equation is obtained based on the first-order scattering elastic wave phase-steady method. Based on the perturbation term in the scattering function, the perturbation stiffness matrix of the HTI medium is obtained; Substituting the perturbation stiffness matrix of the HTI medium into the relationship between the scattering function and the reflection characteristic equation of the HTI medium, an approximate formula for the seismic reflection coefficient of the anisotropic parameters of the HTI medium is obtained. The anisotropy parameters of the HTI medium include: the bulk modulus K of the rock matrix. m Shear modulus μ of the rock matrix m The rock's density ρ, porosity φ, effective pressure P, and fracture normal weakness δ N and crack tangential weakness δ T .

3. The method for predicting pre-stack cracks under velocity anisotropy constraints according to claim 2, characterized in that, The relationship between the scattering function and the HTI medium reflection characteristic equation is as follows: In the formula, R pp (θ) represents the approximate formula for the reflection coefficient of PP waves, ρ b S represents the background medium density. PP (r0) represents the scattering function value at point r = r0, where r represents the position of the scattering function and r0 represents the initial position of the scattering function.

4. The method for predicting pre-stack cracks under velocity anisotropy constraints according to claim 2, characterized in that, The disturbance stiffness matrix of the HTI medium is: in, ΔC 44 =-qΔφμ m +Dm m ΔC 55 =-qΔφμ m -m m Dd T +Dm m , In the formula, ΔC HTI The symbol Δ represents the disturbance stiffness of the HTI medium, indicating the amount of disturbance to the parameter; φ represents the porosity of the rock, Δφ represents the amount of porosity disturbance to the rock; P represents the effective pressure of the rock, ΔP represents the amount of effective pressure disturbance to the rock; K m ΔK represents the bulk modulus of the rock matrix. m δ represents the disturbance of the bulk modulus of the rock matrix. N With δ T Both represent the rock physical parameters of the fracture in the Schoenberg linear slip model, namely the fracture normal weakness and fracture tangential weakness; Δδ N With Δδ T μ represents the weak perturbation in the normal direction and the weak perturbation in the tangential direction of the crack, respectively. m Δμ represents the shear modulus of the rock matrix. m This represents the amount of shear modulus disturbance in the rock matrix; x represents the ratio of the Lamé parameter to the bulk modulus of the rock. M = λ + 2μ, where M represents the bulk modulus of the rock, μ represents the shear modulus of the rock, and λ represents the Lamé parameter.

5. A method for predicting pre-stack cracks under velocity anisotropy constraints according to claim 1 or 2, characterized in that, The approximate formula for the seismic reflection coefficient of the anisotropic parameters of the HTI medium is as follows: in, In the formula, The reflection coefficient of a fractured reservoir varies with the incident angle θ and azimuth angle. The change in θ, where θ represents the angle of incidence, i.e., the angle of incidence of the seismic wave. The azimuth angle represents the angle between the azimuth angle of the observation line and the axis of symmetry of the crack. φ Observe Indicates the observable azimuth angle, φ sym Indicates the symmetrical azimuth angle of the crack; The bulk modulus K of the rock matrix m The coefficient K that varies with the incident angle θ m0 a represents the mean bulk modulus of the rock matrix. μm (θ) represents the shear modulus μ of the rock matrix. m The coefficient a that varies with the incident angle θ ρ (θ) represents the coefficient by which the rock density ρ varies with the incident angle θ, Δρ represents the disturbance of the rock density, ρ0 represents the mean rock density, and a φ (θ) represents the coefficient by which the porosity of the rock varies with the incident angle θ, φ0 represents the mean porosity of the rock, and a p (θ) represents the coefficient of variation of the effective pressure P of the rock with the incident angle θ, and P0 represents the mean effective pressure of the rock. This indicates that the crack normal weakness varies with the incident angle θ and azimuth angle. The coefficient of change This indicates that the tangential weakness of the crack varies with the incident angle θ and the azimuth angle. The coefficient of variation, g1, represents the ratio of the bulk modulus of the rock matrix to the longitudinal wave modulus of the rock, i.e., g1 = K. m / M; g2 represents the ratio of the shear modulus of the rock matrix to the longitudinal wave modulus of the rock, i.e., g2 = μ m / M; g3 represents the ratio of the rock's shear modulus to its longitudinal wave modulus, i.e. Where M and μ represent the longitudinal wave modulus and shear modulus of the rock, respectively.

6. The method for predicting pre-stack cracks under velocity anisotropy constraints according to claim 1, characterized in that, The relationship between the anisotropy parameters and velocity azimuth of the HIT medium is obtained according to the approximate formula for the reflection coefficient. Then, through the inversion of the velocity azimuth anisotropy parameters, a velocity azimuth anisotropy-constrained seismic inversion objective function for the HIT medium is established, including: Based on the aforementioned reflection coefficient approximation formula, the fracture rock physical parameters δ of the Schoenberg linear slip model are obtained. N and δ T ; Using the anisotropic parameters of the HTI medium to be inverted as the parameters to be inverted, the velocity orientation of the HTI medium anisotropic parameters is inverted to obtain the HTI medium anisotropic parameter matrix. Based on the physical parameters of the fractured rock and the anisotropic parameter matrix of the HTI medium using the Schoenberg linear slip model, the velocity orientation and anisotropic constraints are obtained. The objective function for seismic inversion of the HIT medium is established based on the velocity orientation and anisotropy constraints.

7. The method for predicting pre-stack cracks under velocity anisotropy constraints according to claim 6, characterized in that, The physical parameters δ of the fractured rock in the Schoenberg linear slip model are obtained based on the aforementioned approximation formula for the reflection coefficient. N and δ T ,include: Based on the aforementioned approximate formula for the reflection coefficient, the seismic P-wave phase velocity of the HTI medium is characterized using the Thomsen anisotropy parameter. Based on the relationship between the Schoenberg linear slip model and the Hudson thin coin-shaped fracture model, the relationship between the Thomsen anisotropy parameters and the fracture rock physical parameters of the Schoenberg linear slip model was obtained by inversion. Based on the relationship between the Thomsen anisotropy parameters and the fracture rock physical parameters, and the seismic P-wave phase velocity, the fracture rock physical parameters δ of the Schoenberg linear slip model are obtained. N and δ T .

8. The method for predicting pre-stack cracks under velocity anisotropy constraints according to claim 7, characterized in that, The seismic P-wave phase velocity of the HIT medium is: In the formula, Let θ represent the incident angle of the seismic wave and the included angle be θ. The phase velocity of seismic P-waves in HIT media, V P0 θ represents the phase velocity of the seismic P-wave, and θ represents the angle of incidence, i.e., the angle of incidence of the seismic wave. The angle represents the angle between the azimuth of the observation line and the axis of symmetry of the crack; where... φ Observe φ represents the observable azimuth angle. sym Indicates the symmetrical azimuth angle of the crack; V iso ε represents the isotropic velocity of the seismic P-wave. (V) δ (V) and γ (V) All represent the Thomsen anisotropy parameters in the VTI background.

9. A method for predicting pre-stack cracks under velocity anisotropy constraints according to claim 7 or 8, characterized in that, The relationship between the Thomsen anisotropy parameters and the fracture rock physical parameters of the Schoenberg linear slip model is expressed as follows: In the formula, g represents the square of the ratio of seismic transverse to longitudinal wave velocities under isotropic background. V S Indicates the transverse wave velocity of the rock. δ represents the phase velocity of the seismic P-wave in the HIT medium. N and δ T All represent the rock physical parameters of the Schoenberg fracture model; The expression for the fracture rock physical parameters of the Schoenberg linear slip model is as follows: The velocity orientation and anisotropy constraints are as follows: P vel =λ i (or i -l i m i ) Τ (or i -l i m i ), In the formula, λ i δ N and δ T constraint coefficient, in, δ N constraint coefficient, δ T constraint coefficient; η i η represents the scaling factor between the parameters to be inverted and the initial values. i =1 / 2×ln(m) i / m i0 ), m i0 This represents the initial value of the anisotropy parameter matrix of the HIT medium obtained from velocity azimuth inversion; l i This represents the integral operator in a mathematical formula. t i t0 represents the sampling time of the i-th sampling point in the observed seismic record, where i represents the sequence number; t0 represents the initial sampling time.

10. A method for predicting pre-stack cracks under velocity anisotropy constraints according to any one of claims 1-9, characterized in that, The objective function for seismic inversion of the HIT medium is: In the formula, Let ΔS represent the seismic inversion objective function for the HIT medium, ΔG′ represent the perturbation of the seismic record matrix, and ΔG′ represent the forward modeling operator matrix for azimuth amplitude difference. This represents the reflection coefficient matrix of the parameters to be inverted after decorrelation. The variance of random noise is represented by K, the total number of seismic traces involved in the inversion is represented by i, and m represents the index. i This represents the anisotropy parameter matrix of the HIT medium obtained from velocity azimuth inversion. Π represents the variance of the parameter to be inverted. vel This indicates the velocity orientation and anisotropy constraints.

11. The method for predicting pre-stack cracks under velocity anisotropy constraints according to claim 1, characterized in that, The inversion solution of the seismic inversion objective function of the HIT medium to obtain the anisotropic parameters of the target HIT medium characterizing the anisotropic intensity includes: The anisotropy parameters of the HTI medium to be inverted are used as the parameters to be inverted, and the posterior probability density distribution function based on the parameters to be inverted and the observed seismic records is obtained. The objective function for seismic inversion of the HIT medium is solved based on the posterior probability density distribution function to obtain the parameter matrix m to be inverted; Based on the simplified approximate formula for the reflection coefficient, the forward modeling equation of the HTI medium observation seismic record characterized by the anisotropy parameters of the HTI medium is obtained, and the observation seismic record matrix S is obtained according to the forward modeling equation of the HTI medium observation seismic record. The observed seismic record matrix S is transformed using the decorrelated coefficient matrix G′ and the decorrelated parameter matrix m′ to be inverted, to obtain the decorrelated observed seismic record matrix S′. Based on the decorrelation coefficient matrix G′ and the decorrelation observed seismic record matrix S′, the objective function for seismic inversion of the HIT medium is optimized to obtain the objective function for minimum azimuth amplitude difference. The objective function for minimizing azimuth amplitude difference is further solved to obtain the anisotropic parameter δ of the target HTI medium characterizing the anisotropic intensity. N and δ T .

12. The method for predicting pre-stack cracks under velocity anisotropy constraints according to claim 11, characterized in that, The process of obtaining the posterior probability density distribution function based on the parameters to be inverted and the observed seismic records includes: Assuming the parameters to be inverted follow a Cauchy distribution, we obtain the probability distribution of the parameters to be inverted. Assuming that the noise in the observed seismic records follows a Gaussian distribution, determine the likelihood function p(S|m) between the observed seismic records and the parameter matrix to be inverted; Based on the probability distribution of the parameters to be inverted and the likelihood function, the posterior probability density distribution function based on the parameters to be inverted and the observed seismic records is obtained.

13. The method for predicting pre-stack cracks under velocity anisotropy constraints according to claim 12, characterized in that, The probability distribution of the parameters to be inverted is as follows: In the formula, p Cauchy (m) represents the probability distribution of the parameters to be inverted, W represents the number of sampling points in the observed seismic records, and i represents the sequence number; R represents the variance of the parameter vector m to be inverted. i Let R represent the reflection coefficient corresponding to the i-th sampling point, and let R represent the reflection coefficient of the parameter vector m to be inverted. The likelihood function is: In the formula, p(S|m) represents the likelihood function of the observed seismic record and the parameter matrix to be inverted, and σ k This represents the noise distribution in the observed seismic record, where S represents the observed seismic record containing noise, i.e., the observed seismic record. denoted by , where represents the variance of the noise distribution in the observed seismic record, G represents the joint matrix of the wavelet matrix and the coefficient matrix of the parameter to be inverted, and m represents the parameter matrix to be inverted; By combining the probability distribution formula of the parameter to be inverted with the likelihood function formula, the posterior probability density distribution function is obtained.

14. A method for predicting pre-stack cracks under velocity anisotropy constraints according to any one of claims 11-13, characterized in that, The posterior probability density distribution function is: In the formula, W represents the number of sampling points in the observed seismic records.

15. The method for predicting pre-stack cracks under velocity anisotropy constraints according to claim 11, characterized in that, The forward modeling equations for HTI medium observation seismic records, characterized by anisotropic parameters of the HTI medium, are obtained based on the simplified approximate formula for the reflection coefficient, including: The relative difference of the parameters to be inverted in the approximate formula of the reflection coefficient is expressed as the difference of their natural logarithms, so as to remove the fractional terms in the approximate formula of the reflection coefficient and obtain the simplified approximate formula of the reflection coefficient. By convolving the simplified approximate formula for the reflection coefficient with the seismic wavelet, the forward modeling equation for HTI medium seismic records, characterized by the anisotropic parameters of the HTI medium, is obtained.

16. A method for predicting pre-stack cracks under velocity anisotropy constraints according to claim 11 or 15, characterized in that, The simplified approximate formula for the reflection coefficient is: In the formula, R pp (θ,φ) represents the reflection coefficient of a fractured reservoir as a function of the incident angle θ and azimuth angle. Changes, The bulk modulus K of the rock matrix m The coefficient that varies with the incident angle θ. The shear modulus μ of the rock matrix m The coefficient a that varies with the incident angle θ ρ (θ) represents the coefficient that represents the density ρ of the rock as a function of the incident angle θ. φ (θ) represents the coefficient that represents the variation of the rock's porosity φ with the incident angle θ, a p (θ) represents the coefficient that represents the effective pressure P of the rock as a function of the incident angle θ. Indicates the crack normal weakness δ N With incident angle θ and azimuth angle The coefficient of change Indicates the tangential weakness δ of the crack T With incident angle θ and azimuth angle The coefficient of change The bulk modulus K of the rock matrix m The contribution to the reflection coefficient, K m Indicates the bulk modulus of the rock matrix; The shear modulus μ of the rock matrix m Contribution to the reflection coefficient μ m R represents the shear modulus of the rock matrix. ρ R represents the contribution of rock density ρ to the reflectance coefficient. ρ =Δln(ρ), where ρ represents the density of the rock; R φ R represents the contribution of rock porosity φ to the reflection coefficient. φ =Δln(φ), where φ represents the porosity of the rock; R p R represents the contribution of the effective pressure P of the rock to the reflection coefficient. P =Δln(P), where P represents the effective pressure of the rock; Indicates the crack normal weakness δ N The contribution to the reflection coefficient, Indicates the tangential weakness δ of the crack T The contribution to the reflection coefficient, δ N With δ T All represent the rock physical parameters of the fracture in the Schoenberg linear slip model, namely the fracture normal weakness parameter and the fracture tangential weakness parameter; The forward modeling equation for HTI medium seismic records, characterized by the anisotropy parameters of the HTI medium, is as follows: In the formula, S represents the observed seismic record containing noise. This represents seismic wavelets from different directions and angles.

17. The method for predicting pre-stack cracks under velocity anisotropy constraints according to claim 11, characterized in that, The process of obtaining the observed seismic record matrix S based on the forward modeling equation of the HTI medium observed seismic record includes: Through isotropic background earthquake records S iso With anisotropic earthquake records The observed seismic records in the forward modeling equation of the HTI medium are represented in the form of a sum, and the observed seismic record matrix S is obtained.

18. The method for predicting pre-stack fractures under velocity anisotropy constraints according to claim 17, wherein the method uses isotropic background seismic records S iso With anisotropic earthquake records The observed seismic records in the forward modeling equation of the HTI medium are represented in the form of a sum, resulting in the observed seismic record matrix S, which includes: Through isotropic background earthquake records S iso With anisotropic earthquake records The summation form of the observed seismic records in the forward modeling equations of the HTI medium seismic records is expressed as: In the formula, S represents the observed seismic record containing noise, S iso This represents isotropic background earthquake records. Represents anisotropic earthquake records; The observed seismic records, after being represented, are converted into matrix form to obtain the observed seismic record matrix as follows: In the formula, [S] represents the observed seismic record matrix, G represents the joint matrix of the wavelet matrix and the coefficient matrix of the parameters to be inverted, and G = Wa(θ); G iso The joint coefficient matrix of the isotropic background parameters is represented. Let m represent the joint coefficient matrix of the parameters to be inverted. iso This represents the parameter matrix to be inverted against an isotropic background. The bulk modulus K of the rock matrix m Contribution to the reflection coefficient The shear modulus μ of the rock matrix m The contribution of R to the reflection coefficient ρ R represents the contribution of rock density ρ to the reflectance coefficient. φ R represents the contribution of rock porosity φ to the reflection coefficient. p This represents the contribution of the effective pressure P of the rock to the reflection coefficient. This represents the anisotropic parameter matrix to be inverted. Indicates the crack normal weakness δ N The contribution to the reflection coefficient, Indicates the tangential weakness δ of the crack T The contribution to the reflection coefficient.

19. The method for predicting pre-stack cracks under velocity anisotropy constraints according to claim 11, characterized in that, The process of transforming the observed seismic record matrix S using the decorrelated coefficient matrix G′ and the decorrelated parameter matrix m′ to obtain the decorrelated observed seismic record matrix S′ includes: Based on Bayesian theory, the isotropic forward modeling operator G is obtained using the AVAZ inversion method for anisotropic media. iso and anisotropic orthogonal operators Using the covariance matrix C c For the isotropic forward operator G iso Anisotropic orthogonal operators The parameter matrix m to be inverted is then subjected to decorrelation processing to obtain the decorrelation coefficient matrix G′ and the decorrelation parameter matrix m′ to be inverted. The observed seismic record matrix S is transformed using the decorrelation coefficient matrix G′ and the decorrelation parameter matrix m′ to be inverted, and the decorrelation observed seismic record matrix S′ is obtained.

20. The method for predicting pre-stack cracks under velocity anisotropy constraints according to claim 19, characterized in that, The isotropic forward operand G iso for: In the formula, [G iso ] NMF×5F G represents the isotropic forward modeling operator when the azimuth seismic gather has N incident angles, M azimuth angles, and F reflecting interfaces. iso W represents the wavelet matrix. This represents the matrix of bulk modulus coefficients of the rock matrix at N incident angles. A represents the shear modulus coefficient matrix of the rock matrix at N incident angles. ρ (θ) represents the density coefficient matrix of the rock at N incident angles, A φ (θ) represents the porosity coefficient matrix of the rock at N incident angles, A P (θ) represents the effective pressure coefficient matrix of the rock at N incident angles; The anisotropic forward modeling operator for: In the formula, This represents the anisotropic forward modeling operator when the azimuth seismic gather has N incident angles, M azimuth angles, and F reflecting interfaces. This represents the matrix of crack normal weakness coefficients for N incident angles and M azimuth angles. This represents the crack tangential weakness coefficient matrix with N incident angles and M azimuth angles; The decorrelation coefficient matrix G′ and the parameter matrix m′ to be inverted are obtained by the following formula: In the formula, G′ iso G represents the de-correlated isotropic forward modeling operator. iso This represents an isotropic forward modeling operator. This indicates the de-correlated anisotropic forward modeling operator. Let m' represent the anisotropic forward modeling operator. iso m represents the parameter matrix to be inverted from the isotropic background to be removed. iso This represents the parameter matrix to be inverted against an isotropic background. This represents the anisotropic parameter matrix to be inverted to remove correlation. Let represent the anisotropic parameter matrix to be inverted, and u represent the eigenvector.

21. A method for predicting pre-stack cracks under velocity anisotropy constraints according to claim 19 or 20, characterized in that, The decorrelation-free seismic record matrix is ​​as follows: In the formula, [S′] represents the decorrelated observation seismic record matrix.

22. The method for predicting pre-stack cracks under velocity anisotropy constraints according to claim 11, characterized in that, The objective function for minimizing the azimuth amplitude difference is: In the formula, λ represents the minimum azimuth amplitude difference coefficient, ΔG' represents the azimuth amplitude difference forward modeling operator matrix, and λ represents the minimum azimuth amplitude difference coefficient. c Q represents the ratio of the variance of the noise distribution in the observed seismic record to the variance of the parameters to be inverted. c This represents the diagonal matrix of variance coefficients of the inversion parameters. δ N constraint coefficient, δ N The probability distribution function, δ T constraint coefficient, δ T The probability distribution function, ΔS represents the contribution of the parameters to be inverted to the observed seismic record, and ΔS represents the observed seismic record synthesized from the azimuth amplitude difference, i.e. the amount of disturbance to the observed seismic record. in, This represents the variance of the noise distribution in observed seismic records. This represents the variance of the parameters to be inverted, i.e., the variance of the anisotropy parameters of the HTI medium to be inverted; m n This represents the weighting coefficient of the parameter m to be inverted. This represents the variance of the parameter m to be inverted.

23. A pre-stack crack prediction system under velocity anisotropy constraints, characterized in that, include: Establish a unit to obtain an approximate formula for the seismic reflection coefficient of the HTI medium anisotropy parameters; The inversion unit is used to obtain the relationship between the anisotropy parameters of the HTI medium and the velocity azimuth according to the approximate formula of the seismic reflection coefficient, and to establish the objective function of the HIT medium seismic inversion constrained by the velocity azimuth anisotropy through the inversion of velocity azimuth anisotropy. The solution unit is used to perform amplitude difference anisotropy inversion solution on the seismic inversion objective function of the HIT medium, obtain the target HTI medium anisotropy parameters characterizing the anisotropy intensity, and realize pre-stack crack prediction based on velocity azimuth anisotropy constraints.

24. An electronic device, characterized in that, The method includes a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program or instructions to implement the steps of the pre-stack crack prediction method under velocity anisotropy constraints as described in any one of claims 1-21.

25. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program or instructions, which, when executed by a processor, implement the steps of the pre-stack crack prediction method under velocity anisotropy constraints as described in any one of claims 1-22.

26. A computer program product comprising a computer program or instructions, characterized in that, When the computer program or instructions are executed by the processor, they implement the steps of the pre-stack crack prediction method under velocity anisotropy constraints as described in any one of claims 1-22.

Images