Orthogonal fractured reservoir frequency variation AVO inversion and parameter prediction method

Abstract

The invention discloses an orthogonal fractured reservoir frequency variation AVO inversion and parameter prediction method, and relates to the technical field of exploration geophysics, and the method comprises the steps: deducing and reconstructing a fluid saturation horizontal-vertical orthogonal fractured medium rock physical model based on a multi-group fracture-pore medium Chapman model, realizing fluid saturation rock physical modeling containing two groups of orthogonal fracture media; based on an anisotropic reflectivity algorithm and a fluid saturation horizontal-vertical orthogonal fractured medium rock physical model, realizing frequency variation AVO forward modeling of the fractured reservoir, and calculating a frequency variation AVO reflection coefficient and a synthetic seismic record of the fractured reservoir under the interface reflection model; a frequency-varying AVO inversion theory based on a Bayesian framework is constructed by taking a forward modeling process of the frequency-varying AVO as a drive, and the frequency-varying AVO inversion based on the Bayesian inversion framework is realized and reservoir parameter prediction is carried out by utilizing pre-stack seismic AVO data and introducing a simulated annealing-particle swarm algorithm in an inversion process. The method provides a new method and a new means for fine description and prediction of the fractured reservoir.

Description

Technical Field

[0001] The invention relates to the field of exploration geophysical technology, and in particular to an orthogonal fracture type reservoir frequency-variable AVO inversion and parameter prediction method. Background Art

[0002] The dynamic equivalent medium theory model of fluid flow between fractures and pores can more comprehensively describe the frequency-dependent anisotropy characteristics in actual fracture reservoirs. Using the AVO (Amplitude Versus Offset, reflection amplitude varies with shot offset) characteristics of prestack seismic data to identify fracture systems, predict fracture parameters and identify fluid types has become an important means of reservoir fracture prediction. In fractured reservoirs, due to the influence of velocity dispersion and energy attenuation, the elastic parameters of the rock will change with frequency, causing the reflection coefficient to also change with frequency. This phenomenon of reflection coefficient changing with incident angle and frequency is called frequency-dependent AVO characteristics. Based on the dispersive medium rock physics model, the calculation of the frequency-dependent AVO response of fractured reservoirs and the parameter sensitivity analysis constitute the theoretical basis for reservoir parameter prediction using frequency-dependent AVO inversion. At present, frequency-dependent AVO inversion is mainly divided into two categories: data-driven inversion and model-driven inversion. The data-driven frequency-variable AVO inversion theory has been widely applied to shale gas reservoirs, carbonate reservoirs, tight sandstone reservoirs and natural gas hydrate reservoirs, and has achieved good results. However, the theoretical basis of this type of method is the AVO approximate formula, which only considers the case of single interface reflection, and does not consider the layered structure and complex lithological information of the formation. The model-driven frequency-variable AVO inversion can consider one or more dispersion mechanisms and complex formation structures in fractured reservoirs. Existing studies have shown that the model-driven frequency-variable AVO inversion method has significant advantages in predicting complex fractured reservoirs and reservoirs containing multiple dispersion mechanisms, but the above studies only perform inversion estimation on a single set of fracture models.

[0003] Therefore, it is urgent to study a method to quantitatively describe the fracture parameters of multi-group fracture-type reservoirs containing fluids. Summary of the invention

[0004] In order to solve the above technical problems, the present invention discloses a frequency-variable AVO inversion and parameter prediction method for orthogonal fracture reservoirs, which aims to quantitatively describe the parameters of fluid-saturated orthogonal fracture reservoirs, construct a fluid-saturated orthogonal fracture rock physics model and a frequency-variable AVO forward modeling process, provide theoretical support for the technical method of reservoir description and inversion based on dispersion attributes, and establish a frequency-variable AVO inversion process based on the Bayesian inversion framework to achieve frequency-variable AVO inversion prediction of complex fluid and fracture systems. This method can not only improve the accuracy of predicting reservoir parameters, but also quantify uncertainty, provide more comprehensive information for reservoir evaluation, and provide new methods and means for the detailed description and prediction of fracture reservoirs.

[0005] To achieve the above object, the present invention adopts the following technical solutions:

[0006] A frequency-variable AVO inversion and parameter prediction method for orthogonal fracture reservoirs includes the following steps:

[0007] S1. Based on the Chapman model of multiple groups of fracture-pore media, the rock physics model of fluid-saturated horizontal-vertical orthogonal fracture media is derived and reconstructed to achieve rock physics modeling of fluid-saturated media containing two groups of orthogonal fractures;

[0008] S2. Based on the anisotropic reflectivity algorithm and the rock physics model of fluid-saturated horizontal-vertical orthogonal fracture media, the frequency-variant AVO forward modeling of fractured reservoirs is realized, and the frequency-variant AVO reflection coefficient and synthetic seismic records of fractured reservoirs under the interface reflection model are calculated;

[0009] S3. Driven by the forward modeling process of frequency-variable AVO, a frequency-variable AVO inversion theory based on the Bayesian framework is constructed. Using prestack seismic AVO data, the simulated annealing-particle swarm algorithm (SA-PSO) is introduced into the inversion process to realize frequency-variable AVO inversion based on the Bayesian inversion framework and predict reservoir parameters.

[0010] In the model of Chapman (2009), two fracture scales are considered: microscopic spherical pores and mesoscopic fractures. The model of the present invention selects two groups of arbitrarily distributed fractures to be defined as mutually orthogonal and the fracture fluid is a saturated model of oil. In step S1, the specific steps are as follows:

[0011] According to Chapman 92009), a rock physics model describing the velocity dispersion and attenuation characteristics in multiple fracture media is given. The model assumes that the pore radius is equal to the particle size, the fracture size is much larger than the particle size but smaller than the wavelength of the seismic wave, the rock skeleton is isotropic, and the interaction between fractures is ignored. Assuming that there are two groups of orthogonal fractures f1 and f2 in the rock, the polarization angle θ is used. m and azimuth Two angles characterize the distribution of cracks, m = 1, 2 represents two groups of cracks, [X1, X2] crack plane is the horizontal crack plane, [X1, X3] is the vertical crack plane, and it is defined that crack f1 is on the horizontal crack plane and crack f2 is on the vertical crack plane, that is, θ1 = 0°, θ2=90°, Then the normal direction of the crack surface is defined as:

[0012] n 1 =(0,0,1) (1);

[0013] n 2 =(1,0,0) (2);

[0014] The frequency-dependent effective stiffness matrix of the porous equivalent medium with orthogonal cracks is derived as follows:

[0015]

[0016] In the formula, C ij (ω) is the equivalent elastic tensor of rock, is the elastic tensor of the isotropic background matrix, is the elastic correction term caused by porosity, is the elastic correction term caused by horizontal cracks, is the elastic correction term caused by vertical cracks;

[0017] Elastic tensor (stiffness coefficient matrix) of the isotropic background matrix Expressed by elastic parameters Lame coefficients λ and μ, the expression is:

[0018]

[0019] Elastic correction due to porosity The standard matrix form of is:

[0020]

[0021] In the formula, the expressions of d and e are:

[0022]

[0023] In the formula, is the porosity, u is the Poisson's ratio, is the stress tensor of the pore, p* is the fluid pressure in the pore;

[0024] Pore ​​stress tensor The expressions of the fluid pressure p* in the pores are:

[0025]

[0026] In the formula, is the porosity of group m fractures, τ m is the relaxation time scale parameter of m groups of cracks, is the critical pressure of m groups of cracks, K p is the bulk modulus of the pore, is the bulk modulus of m groups of cracks;

[0027] Bulk modulus K of pores p , the bulk modulus of m groups of cracks and the critical pressure of group m fractures The expressions are:

[0028]

[0029] In the formula, κ f is the bulk modulus of the fluid;

[0030] Elastic correction due to horizontal cracks and the elastic correction term due to vertical cracks The elastic parameter matrix of is characterized by the fluid pressure in the fracture and is expressed as:

[0031]

[0032]

[0033] In the formula, the matrix elements The expression is:

[0034]

[0035] In the formula, r m is the aspect ratio of the crack, is the fluid pressure in group m fractures;

[0036] m=1,θ1=0°, For the first group of fractures, the fluid pressure of the fractures arranged horizontally is expressed as:

[0037]

[0038] m=2,θ2=90°, For the second group of fractures, the fluid pressure of the fractures arranged vertically is expressed as:

[0039]

[0040] According to formula (3), the vertical P-wave velocity V of porous rock with orthogonal fractures saturated with fluid is obtained as follows: p (ω), vertical shear wave velocity V s(ω) and the attenuation factor The expression is:

[0041]

[0042] In the formula, Re(C 33 (ω)) and Im(C 33 (ω)) represent the real part and imaginary part of the complex number respectively.

[0043] Optionally, in step S2, the saturated fluid orthogonal fracture porous rock physical model is derived into a frequency-dependent anisotropic elastic model, and the equation is simplified by combining the anisotropic reflectivity method. After Fourier transformation, the expression of the anisotropic reflectivity method based on the differential equation system is obtained as follows:

[0044]

[0045] Where z is the depth; b is the displacement-stress vector; ω is the angular frequency; the matrix A is a 6×6 system matrix, which is composed of the elastic tensor C of the corresponding layer ij (ω), density ρ and horizontal slowness components in the x and y directions indicate that, due to the different elastic parameters of different media, the eigenvalues ​​and eigenvectors of the system matrix A will vary with the media;

[0046] To solve formula (31), based on the condition that displacement and stress are continuous on any plane, the anisotropic reflectivity method introduces the propagation matrix operator and calculates the reflection and transmission coefficient matrix through recursive equations. Based on the saturated fluid orthogonal fracture porous rock physics model, the upgoing wave reflection coefficient matrix R is re-derived using the eigenvalues ​​and eigenvectors of the system matrix A. U , downlink wave reflection coefficient matrix R D , upgoing wave transmission coefficient matrix T U and the downlink transmission coefficient matrix T D , the expression is:

[0047]

[0048] Where, the matrix elements are the frequency domain reflection responses; For example, the first subscript in the element is the type of reflected wave, the second subscript is the type of incident wave, and the p, 1, and 2 in the subscript represent q p Wave, q s1 Wave and Q s2 Waves Three types of waves;

[0049] When the energy is incident from top to bottom, the reflection and transmission coefficient matrix R of the anisotropic medium is expressed as:

[0050]

[0051] In the formula, R U , R D 、T U 、T D are the reflection coefficient matrix and transmission coefficient matrix of the upgoing and downgoing waves;

[0052] The propagation matrix operator is introduced into the anisotropic reflectivity method and the reflection and transmission coefficient matrices are calculated through recursive equations. The reflection and transmission coefficient matrices of the i+1th layer are calculated by the i-th layer recursive relationship. The recursive relationship equation is:

[0053]

[0054] In the formula, r U 、r D ,t U ,t D The interface depth z=z i+1 The reflection matrix and transmission matrix at E U is the eigenvector of the upgoing wave reflection coefficient, E D is the eigenvector of the downlink wave reflection coefficient;

[0055] E U 、E D The eigenvectors are:

[0056]

[0057] In the formula, h = z i+1 -z is the depth difference between the i-th layer and the i+1-th layer; λ is the eigenvalue of the corresponding system matrix A

[0058] The frequency-variable AVO reflection coefficient matrix in the entire layered medium is obtained by iterative calculation using equations (37)-(40); the reflection matrix includes the effects of incident P-wave reflection, transmission, converted waves, and multiple waves. The time domain prestack synthetic seismic record is calculated by the inverse Fourier transform of the frequency-variable AVO reflection coefficient and the seismic wavelet spectrum, and the formula is:

[0059]

[0060] In the formula, S PP is the pre-stack synthetic seismic record in the time domain; R PP (ω, θ) is the frequency-varying AVO reflection coefficient matrix, θ is the incident angle of the seismic wavelet, and F(ω) is the seismic wavelet spectrum.

[0061] Optionally, in step S3, the core of the Bayesian inversion theory is based on the idea of ​​probability statistics, that is, using prior information and likelihood function to obtain posterior information, so that the problem of obtaining the minimum value of the inversion objective function is transformed into the problem of maximizing the posterior probability distribution of the model parameters. In order to realize the prediction of reservoir parameters, based on the frequency-variable AVO forward modeling process, the posterior probability formula of prestack seismic data for reservoir parameters under the Bayesian inversion framework is derived as follows:

[0062]

[0063] In the formula, m(ε f1 , ε f2 , Φ p , a f1 , a f2 ) is the model reservoir parameter; d is the model observation data (prestack seismic data); P(m)) is the prior information of the model reservoir parameter; P(d|m)) is the likelihood function; ∫P(d|m))·P(m))dm is the probability normalization factor, which is a constant and is used to normalize the final inversion result;

[0064] The observed data d and the forward synthetic seismic record S are calculated by the calculation model PP The error between them, the objective function is constructed as:

[0065] ΔE=∑||d-SP P ||2 (45);

[0066] By exponentially transforming the error function, the likelihood function P(d|m(ε f1 , ε f2 , Φ p , a f1 , a f2 )):

[0067] P(d|m(ε f1 , ε f2 , Φ p , a f1 , a f2 ))∝a·exp(-b,ΔE) (46);

[0068] Where b is a constant;

[0069] Obtaining prior information P(m(ε) from logging information f1 , ε f2 , Φ p , a f1 , a f2 ), using the model observation data as constraints, the reservoir parameter m(ε f1 , ε f2 , Φ p, a f1 , a f2 ), the objective function is constructed as:

[0070] m(ε f1 ,ε f2 ,Φ p ,a f1 ,a f2 )=arg max[P(m(ε f1 ,ε f2 ,Φ p ,a f1 ,a f2 )|d)](47);

[0071] The simulated annealing particle swarm algorithm is used to solve the objective function formula (47). By searching for the parameter m(ε f1 , ε f2 , Φ p , a f1 , a f2 ) is the optimal solution, so that the posterior probability P(m(ε f1 , ε f2 , Φ p , a f1 , a f2 )|d)maximum; end the calculation and output the inversion results.

[0072] The beneficial effect of the present invention is that the present invention proposes a frequency-variable AVO inversion and parameter prediction method for orthogonal fracture-type reservoirs, performs rock physical description on the dispersion characteristics of fracture-type reservoirs, establishes a frequency-variable AVO forward modeling process for fracture-type reservoirs, and constructs an inversion framework for fine description and prediction of fracture-type reservoirs. First, based on the Chapman model of multiple groups of fracture-pore media, the relationship between the complex frequency-variable elastic stiffness coefficient and the fracture parameters is derived, and the rock physical modeling of fluid-saturated horizontal-vertical orthogonal fracture media is realized, and the velocity dispersion and attenuation characteristics of the model medium are studied. Secondly, based on the anisotropic reflectivity algorithm and the rock physical model of orthogonal fracture media, the frequency-variable AVO forward modeling of fracture-type reservoirs is realized, the frequency-variable AVO reflection coefficient and synthetic seismic records of fracture-type reservoirs are calculated, and the frequency-variable AVO characteristics of the reservoir are studied. Finally, based on the Bayesian inversion framework and driven by the forward modeling process of frequency-varying AVO, the inversion objective function was reconstructed, and the simulated annealing particle swarm algorithm (SA-PSO) was introduced into the inversion process to realize the frequency-varying AVO inversion based on the Bayesian inversion framework. At the same time, the prestack seismic azimuthal AVO data were used to predict the fracture parameters, providing new methods and means for the detailed description and prediction of fractured reservoirs. BRIEF DESCRIPTION OF THE DRAWINGS

[0073] Figure 1 It is a schematic diagram of the rock physics model of horizontal-vertical orthogonal fracture media;

[0074] Figure 2 Define schematic diagrams for reflection and transmission coefficients;

[0075] Figure 3 This is the flow chart of frequency-variant AVO inversion for orthogonal fracture reservoirs based on the Bayesian framework;

[0076] Figure 4 It is the vertical wave velocity dispersion and attenuation factor curve of the horizontal-vertical orthogonal fracture medium model;

[0077] Figure 5 The vertical P-wave velocity dispersion and attenuation factor curves of the orthogonal fracture medium model under different fracture parameter conditions;

[0078] Figure 6 It is a schematic diagram of the single interface reflection model;

[0079] Figure 7 The variation characteristics of PP wave reflection coefficient of single interface reservoir model with incident angle and frequency under different fracture parameter conditions;

[0080] Figure 8 The frequency-dependent seismic AVO response of the single-interface reservoir model under different fracture parameters;

[0081] Fig. 9 Synthetic seismic records for reservoir models;

[0082] Fig.10 is the fracture density of the reservoir model (ε f1 , ε f2 )’s likelihood function;

[0083] Fig.11 are two sets of orthogonal crack densities (ε f1 , ε f2 )’s prior distribution;

[0084] Fig.12 is the fracture density of the reservoir model (ε f1 , ε f2 ), where the true value is represented by a white dot;

[0085] Fig.13 is the posterior distribution of the reservoir model fracture density and porosity inversion under different porosity conditions, where the true value is represented by a white dot;

[0086] Fig.14 Figure 2 is the posterior distribution of the inversion of the fracture density and fracture radius of the reservoir model under different fracture radius conditions, where the true value is represented by white dots. DETAILED DESCRIPTION

[0087] In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, rather than all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention. Therefore, the following detailed description of the embodiments of the present invention provided in the drawings is not intended to limit the scope of the invention claimed for protection, but merely represents the selected embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention.

[0088] A frequency-variable AVO inversion and parameter prediction method for orthogonal fracture reservoirs, such as Figure 3 As shown, the following steps are included:

[0089] S1. Based on the Chapman model of multiple groups of fracture-pore media, the rock physics model of fluid-saturated horizontal-vertical orthogonal fracture media is derived and reconstructed to realize the rock physics modeling of fluid-saturated media containing two groups of orthogonal fractures.

[0090] In the model of Chapman (2009), two fracture scales are considered: microscopic spherical pores and mesoscopic fractures. The model of the present invention selects two groups of arbitrarily distributed fractures to be defined as mutually orthogonal and the fracture fluid is a saturated model of oil. The specific steps of step S1 are as follows:

[0091] According to Chapman (2009), a rock physics model describing the velocity dispersion and attenuation characteristics in a multi-group fracture medium is given. The model assumes that the pore radius is equal to the particle size, the fracture size is much larger than the particle size but smaller than the wavelength of the seismic wave, the rock skeleton is isotropic, and the interaction between fractures is ignored. Assuming that there are two groups of orthogonal fractures f1 and f2 in the rock, the polarization angle θ is used. m and azimuth Two angles characterize the distribution of cracks, m = 1, 2 represents two groups of cracks, [X1, X2] crack plane is the horizontal crack plane, [X1, X3] is the vertical crack plane, and it is defined that crack f1 is on the horizontal crack plane and crack f2 is on the vertical crack plane, that is, θ1 = 0°, θ2=90°, like Figure 1 As shown, the normal direction of the crack surface is defined as:

[0092] n 1 =(0,0,1) (1);

[0093] n 2=(1,0,0) (2);

[0094] The frequency-dependent effective stiffness matrix of the porous equivalent medium with orthogonal cracks is derived as follows:

[0095]

[0096] In the formula, C ij (ω) is the equivalent elastic tensor of rock, is the elastic tensor of the isotropic background matrix, is the elastic correction term caused by porosity, is the elastic correction term caused by horizontal cracks, is the elastic correction term caused by vertical cracks;

[0097] Elastic tensor (stiffness coefficient matrix) of the isotropic background matrix Expressed by elastic parameters Lame coefficients λ and μ, the expression is:

[0098]

[0099] Elastic correction due to porosity The standard matrix form of is:

[0100]

[0101] In the formula, the expressions of d and e are:

[0102]

[0103] In the formula, is the porosity, v is the Poisson's ratio, is the stress tensor of the pore, p* is the fluid pressure in the pore;

[0104] Pore ​​stress tensor The expressions of the fluid pressure p* in the pores are:

[0105]

[0106] In the formula, is the porosity of group m fractures, τ m is the relaxation time scale parameter of m groups of cracks, is the critical pressure of m groups of cracks, K p is the bulk modulus of the pore, is the bulk modulus of m groups of cracks;

[0107] Bulk modulus K of pores p , the bulk modulus of m groups of cracks and the critical pressure of group m fractures The expressions are:

[0108]

[0109] In the formula, κ f is the bulk modulus of the fluid;

[0110] Elastic correction due to horizontal cracks and the elastic correction term due to vertical cracks The elastic parameter matrix of is characterized by the fluid pressure in the fracture and is expressed as:

[0111]

[0112] In the formula, the matrix elements The expression is:

[0113]

[0114]

[0115] In the formula, r m is the aspect ratio of the crack, is the fluid pressure in group m fractures;

[0116] m=1,θ1=0°, For the first group of fractures, the fluid pressure of the fractures arranged horizontally is expressed as:

[0117]

[0118] m=2,θ2=90°, For the second group of fractures, the fluid pressure of the fractures arranged vertically is expressed as:

[0119]

[0120] According to formula (3), the vertical P-wave velocity V of porous rock with orthogonal fractures saturated with fluid is obtained as follows: p (ω), vertical shear wave velocity V s (ω) and the attenuation factor The expression is:

[0121]

[0122] In the formula, Re(C 33 (ω)) and Im(C 33 (ω)) represent the real part and imaginary part of the complex number respectively.

[0123] S2. Based on the anisotropic reflectivity algorithm and the rock physics model of fluid-saturated horizontal-vertical orthogonal fracture media, the frequency-variant AVO forward modeling of fractured reservoirs is realized, and the frequency-variant AVO reflection coefficient and synthetic seismic records of fractured reservoirs under the interface reflection model are calculated.

[0124] When seismic waves propagate in a complex fractured medium containing fluid, it will cause velocity dispersion and energy attenuation within the seismic frequency band, resulting in the reflection coefficient and transmission coefficient of the seismic waves also varying with frequency. The horizontal-vertical orthogonal fracture medium rock physics model constructed by the present invention can derive a frequency-dependent anisotropic elastic model, and the frequency-dependent reflection coefficient of the reservoir can be calculated based on the model and the anisotropic reflectivity method. Compared with the traditional Zoeppritz equation and its approximate solution, the anisotropic reflectivity method solves the problem that it is based on the assumption of vertical incidence or only considers the case of single interface reflection.

[0125] The specific steps include:

[0126] The saturated fluid orthogonal fracture porous rock physics model is derived into a frequency-dependent anisotropic elastic model, and the equation is simplified by combining the anisotropic reflectivity method. After Fourier transformation, the expression of the anisotropic reflectivity method based on the differential equation system is obtained as follows:

[0127]

[0128] Where z is the depth; b is the displacement-stress vector; ω is the angular frequency; the matrix A is a 6×6 system matrix, which is composed of the elastic tensor C of the corresponding layer ij (ω), density ρ and horizontal slowness components in the x and y directions indicate that, due to the different elastic parameters of different media, the eigenvalues ​​and eigenvectors of the system matrix A will vary with the media;

[0129] To solve formula (31), based on the condition that displacement and stress are continuous on any plane, the anisotropic reflectivity method introduces the propagation matrix operator and calculates the reflection and transmission coefficient matrix through recursive equations. Based on the saturated fluid orthogonal fracture porous rock physics model, the upgoing wave reflection coefficient matrix R is re-derived using the eigenvalues ​​and eigenvectors of the system matrix A. U , downlink wave reflection coefficient matrix R D , upgoing wave transmission coefficient matrix T U and the downlink transmission coefficient matrix T D , the expression is:

[0130]

[0131] Where, the matrix elements are the frequency domain reflection responses; For example, the first subscript in the element is the type of reflected wave, the second subscript is the type of incident wave, and the p, 1, and 2 in the subscript represent the three types of waves (q p Wave, q s1 Wave and Q s2 Wave);

[0132] When energy is incident from top to bottom, such as Figure 2 As shown, the reflection and transmission coefficient matrix R of anisotropic media is expressed as:

[0133]

[0134] In the formula, R U , R D 、T U 、T D are the reflection coefficient matrix and transmission coefficient matrix of the upgoing and downgoing waves;

[0135] The propagation matrix operator is introduced into the anisotropic reflectivity method and the reflection and transmission coefficient matrices are calculated through recursive equations. The reflection and transmission coefficient matrices of the i+1th layer are calculated by the i-th layer recursive relationship. The recursive relationship equation is:

[0136]

[0137] In the formula, r U 、r D ,t U ,t D The interface depth z=z i+1 The reflection matrix and transmission matrix at E U is the eigenvector of the upgoing wave reflection coefficient, E D is the eigenvector of the downlink wave reflection coefficient;

[0138] E U 、E D The eigenvectors are:

[0139]

[0140] In the formula, h = z i+1 -z is the depth difference between the i-th layer and the i+1-th layer; λ is the eigenvalue of the corresponding system matrix A

[0141] The frequency-variable AVO reflection coefficient matrix in the entire layered medium is obtained by iterative calculation using equations (37)-(40); the reflection matrix includes the effects of incident P-wave reflection, transmission, converted waves, and multiple waves. The time domain prestack synthetic seismic record is calculated by the inverse Fourier transform of the frequency-variable AVO reflection coefficient and the seismic wavelet spectrum, and the formula is:

[0142]

[0143] In the formula, S PP is the pre-stack synthetic seismic record in the time domain; R PP (ω, θ) is the frequency-varying AVO reflection coefficient matrix, θ is the incident angle of the seismic wavelet, and F(ω) is the seismic wavelet spectrum.

[0144] S3. Driven by the forward modeling process of frequency-variable AVO, a frequency-variable AVO inversion theory based on the Bayesian framework is constructed. Using prestack seismic AVO data, a simulated annealing particle swarm optimization (SA-PSO) algorithm is introduced into the inversion process to realize frequency-variable AVO inversion based on the Bayesian inversion framework and predict reservoir parameters.

[0145] Based on prestack seismic data, the present invention constructs a Bayesian inversion framework for jointly estimating multiple groups of fracture parameters in a fracture reservoir model, wherein all other parameters in the model are considered to be known. Based on the Bayesian inversion framework, driven by the forward modeling process of frequency-variant AVO, the inversion objective function is reconstructed, and the simulated annealing-particle swarm algorithm (SA-PSO) is introduced into the inversion process to realize frequency-variant AVO inversion based on the Bayesian inversion framework and perform reservoir parameter prediction.

[0146] The core of Bayesian inversion theory is based on the idea of ​​probability statistics, that is, using prior information and likelihood function to obtain posterior information, so that the problem of finding the minimum value of the inversion objective function is transformed into the problem of maximizing the posterior probability distribution of the model parameters. In order to realize the prediction of reservoir parameters, based on the frequency-variable AVO forward modeling process, the posterior probability formula of prestack seismic data for reservoir parameters under the Bayesian inversion framework is derived as follows:

[0147]

[0148] In the formula, m(ε f1 , ε f2 , Φ p , a f1 , a f2 ) is the model reservoir parameter; d is the model observation data, prestack seismic data; P(m)) is the prior information of the model reservoir parameter; P(d|m)) is the likelihood function; ∫P(d|m))·P(m))dm is the probability normalization factor, which is a constant and is used to normalize the final inversion result;

[0149] The observed data d and the forward synthetic seismic record S are calculated by the calculation model PP The error between them, the objective function is constructed as:

[0150] ΔE=∑||dS PP ||2 (45);

[0151] By exponentially transforming the error function, the likelihood function P(d|m(ε f1 , ε f2 , Φ p , a f1 , a f2 )):

[0152] P(d|m(ε f1 , ε f2 , Φ p , a f1 , a f2 ))∝a·exp(-b,ΔE)(46);

[0153] Where b is a constant;

[0154] Obtaining prior information P(m(ε) from logging information f1 , ε f2 , Φ p , a f1 , a f2 ), using the model observation data as constraints, the reservoir parameter m(ε f1 , ε f2 , Φ p , a f1 , a f2 ), the objective function is constructed as:

[0155] m(ε f1 ,ε f2 ,Φ p ,a f1 ,a f2 )=arg max[P(m(ε f1 ,ε f2 ,Φ p ,a f1 ,a f2 )|d)](47);

[0156] The simulated annealing particle swarm algorithm is used to solve the objective function formula (47). By searching for the parameter m(ε f1 , ε f2 , Φ p , a f1 , a f2 ) is the optimal solution, so that the posterior probability P(m(ε f1 , ε f2 , Φ p , a f1 , a f2 )|d)maximum;

[0157] End the calculation and output the inversion results.

[0158] Theoretical model analysis

[0159] 1. Analysis of the effects of fracture parameters on velocity dispersion and attenuation characteristics of fluid-saturated horizontal-vertical orthogonal fracture media

[0160] Figure 4 (a) is the vertical wave velocity dispersion curve of the horizontal-vertical orthogonal crack medium model calculated based on the multi-scale Chapman dynamic equivalent medium model. Figure 4 (b) is the attenuation factor curve. The specific model parameters are shown in Table 1. Figure 4 (a) and (b) show that within the seismic frequency band, the vertical P-wave velocity Vp exhibits obvious dispersion and attenuation characteristics, while the vertical S-wave velocity Vs is affected by the frequency-independent parameter C. 55 The influence of does not show a trend of changing with frequency.

[0161] Table 1 Parameters of the theoretical model of horizontal-vertical orthogonal cracks

[0162]

[0163] In order to study the influence of crack parameters on the dispersion characteristics of the model, the present invention sets six different crack parameter models, and other background parameters remain unchanged. By calculating and analyzing the changes in vertical P-wave velocity and quality factor under different crack parameter conditions, the results are as follows: Figure 5 As shown, (a) is the vertical P-wave velocity dispersion curve of the orthogonal fracture medium model under different porosity conditions; (b) is the vertical P-wave velocity attenuation factor curve of the orthogonal fracture medium model under different porosity conditions; (c) is the vertical P-wave velocity dispersion curve of the orthogonal fracture medium model under different fracture density conditions; (d) is the vertical P-wave velocity attenuation factor curve of the orthogonal fracture medium model under different porosity conditions; (e) is the vertical P-wave velocity dispersion curve of the orthogonal fracture medium model under different fracture radius conditions; (f) is the vertical P-wave velocity attenuation factor curve of the orthogonal fracture medium model under different fracture radius conditions. Figure 5 (a), (c), and (e) show that when the porosity, fracture density, and fracture radius of the model change respectively, the dispersion phenomenon of the model is significant. The vertical P-wave velocity first increases with the increase of frequency, and then gradually flattens out and reaches a maximum value when the frequency reaches 100 Hz. Figure 5 (b), (d), and (f) show that the frequency band position of the model attenuation peak moves with the change of fracture parameters. Since there are obvious differences in the dispersion and attenuation of the P-wave velocity of orthogonal fractures under different fracture parameter conditions, the seismic dispersion attributes corresponding to these differences can be used to predict reservoir fracture parameters.

[0164] 2. Frequency-dependent AVO response characteristics of fracture parameters to orthogonal fracture reservoirs

[0165] Before verifying the effectiveness of the inversion method, it is necessary to study the effect of the changes in the fracture parameters to be inverted in the orthogonal fracture reservoir on the PP wave frequency-dependent reflection coefficient and analyze the sensitivity of the PP wave frequency-dependent AVO response. Figure 5 The P-wave velocity dispersion curves with different fracture parameters are shown in the figure. The present invention constructs a single-interface reservoir model with an upper layer of isotropic shale and a lower layer of horizontal-vertical orthogonal fracture dispersion sandstone, such as Figure 6 As shown, the elastic parameters of the shale layer are: Vp = 4.3 km / s, Vs = 2.35 km / s, ρ = 2.68 g / cm 3 , ε=0.2, δ=0.15, γ=0.1. Furthermore, the anisotropic reflectivity method is used to simulate the single interface reservoir model. Assuming that the incident seismic wave has a main frequency of 25 Hz and the azimuth angle is set to 30°, the incident angle range is 0°-30°, and the PP wave reflection coefficient is calculated according to formula (43) with the change of frequency and incident angle and the seismic synthetic record. The results are as follows: Figure 7 and Figure 8 shown. Figure 7 In the figure, (a) is the porosity Variation characteristics of PP wave reflection coefficient of single interface reservoir model with incident angle and frequency under the condition; (b) is the porosity The variation characteristics of the PP wave reflection coefficient of the single interface reservoir model with the incident angle and frequency under the condition; (c) is the fracture density ε f1 = 0.01 under the condition of single interface reservoir model PP wave reflection coefficient with incident angle and frequency; (d) is the fracture density ε f2 = 0.08 under the condition of single interface reservoir model PP wave reflection coefficient with incident angle and frequency; (e) is the fracture radius a f1 =0.05m single interface reservoir model PP wave reflection coefficient with incident angle and frequency; (f) is the fracture radius a f2 =1m single interface reservoir model PP wave reflection coefficient variation characteristics with incident angle and frequency; Figure 8 In the figure, (a) is the porosity Frequency-dependent seismic AVO response of the single-interface reservoir model under the condition of (a) porosity Frequency-dependent seismic AVO response of the single-interface reservoir model under the condition; (c) is the fracture density ε f1 =0.01 under the condition of single interface reservoir model frequency-dependent seismic AVO response; (d) is the fracture density ε f2 = 0.08 single interface reservoir model frequency-dependent seismic AVO response; (e) is the fracture radius a f1 =0.05m single interface reservoir model frequency-dependent seismic AVO response; (f) is the fracture radius a f2=1m single interface reservoir model frequency-dependent seismic AVO response. The results show that in the single interface reservoir model, due to the velocity dispersion and attenuation characteristics corresponding to different porosities, fracture densities and fracture radii (such as Figure 5 shown) are different, resulting in Figure 7 The characteristics of the frequency-dependent reflection coefficient in (a) to (f) are also different. Under different fracture parameters, the PP wave reflection coefficient has a significant frequency-dependent response within the seismic frequency band, and the amplitude of the reflection coefficient increases with the increase of the incident angle. Figure 8 is the corresponding seismic AVO response. It can be seen that when the fracture parameters change, the reflection amplitude also increases with the increase of the incident angle. Figure 7 and Figure 8 It can be considered that the scheme of establishing a more sophisticated fracture reservoir model, calculating the frequency-varying AVO response, and inverting multiple groups of reservoir fracture parameters is feasible.

[0166] 3. Frequency-dependent AVO inversion and parameter prediction of orthogonal fracture reservoirs based on Bayesian framework

[0167] According to the above sensitivity analysis of reservoir fracture parameters, the horizontal-vertical orthogonal reservoir fracture densities are selected as follows: f1 = 0.05 and ε f2 = 0.02, other background parameters remain unchanged, and the forward calculation of single interface reservoir is performed. To simulate the actual earthquake situation, 10% Gaussian noise is added to the synthetic seismic record, and the calculation results are Fig. 9 Synthetic seismic records without noise (a) and with 10% Gaussian noise (b) are then quantitatively inverted using a Bayesian inversion procedure for two sets of fracture density values. Fig.10 For scanning ε f1 and ε f2 The mismatch between the observed data and the model response is calculated by combining the observed data and the model response, thereby calculating the likelihood function P(d|ε f1 , ε f2 ), where (a) is the condition without noise, and (b) is the condition with 10% Gaussian noise. Fig.11 is the corresponding prior probability distribution P(ε f1 , ε f2 ), assuming that the prior information ε f2 is a uniform distribution, ε f1 is a normal distribution. Finally, the posterior probability P(ε f1 , ε f2 |d). Fig.12 The density of two groups of fractures (ε f1 ,ε f2), the results show that in the absence of noise, the two groups of fracture densities obtained by inversion are ε f1 is about 0.05 and ε f2 is about 0.02; in the case of noise, the two groups of fracture densities obtained by inversion are also ε f1 is about 0.05 and ε f2 The inversion results show that despite the presence of noise in the data, the density of the two groups of fractures can still be accurately estimated.

[0168] Furthermore, by changing the prior information of the orthogonal fracture reservoir, the Bayesian inversion steps are repeated to predict the reservoir porosity and fracture radius, thus verifying the effectiveness and applicability of the proposed method. Fig.13 (a) and (b) show the posterior probability results of Bayesian inversion when the reservoir porosity is 10% and 25% respectively, under the condition that the fracture density and porosity are known as prior information. It can be seen from the figure that for reservoirs with larger porosity, the inversion results are more accurate. Fig.14 (a) and (b) show the results of the analysis of the reservoir fracture radius a when the fracture density and fracture radius are known as prior information. f1 =0.05m and a f1 =1m for Bayesian inversion. By comparing the inversion results, it can be found that when the reservoir contains smaller-scale fractures, the inversion method can accurately identify them, but the recognition accuracy of large-scale fractures is weak.

[0169] For fluid-saturated orthogonal fracture reservoirs, the frequency-variable AVO inversion and parameter prediction method for orthogonal fracture reservoirs based on the Bayesian framework proposed in this paper obtains the following conclusions:

[0170] 1) Based on the frequency-dependent dynamic equivalent model with different directions and scales given by the Chapman model, the dispersion and attenuation characteristics of the fracture medium are analyzed when the directions of the two groups of fractures are horizontally and vertically orthogonal. It is found that the changes in fracture parameters (porosity, fracture density and fracture radius) are highly sensitive to velocity dispersion.

[0171] 2) Driven by the equivalent model of fluid-saturated orthogonal fracture reservoir, the frequency-dependent AVO response characteristics under different fracture reservoir parameters were simulated using the anisotropic reflectivity method. By analyzing the changes in the PP wave frequency-dependent reflection coefficient caused by changes in fracture reservoir parameters (such as porosity, fracture density and fracture radius) with frequency and incident angle, the variation law of seismic AVO response was studied. The results show that different fracture parameters have a significant effect on the frequency-dependent AVO response, especially at specific frequencies and incident angles, the change in reflection coefficient is particularly obvious.

[0172] 3) By studying the influence of reservoir fracture parameters on the frequency-variant reflection coefficient of PP waves, a frequency-variant AVO inversion theory based on the Bayesian framework inversion framework driven by the frequency-variant AVO forward modeling process was established. The simulated annealing particle swarm algorithm (SA-PSO) was introduced into the inversion process using prestack seismic AVO data to achieve the inversion of fracture reservoir parameters based on the Bayesian inversion framework. At the same time, the effectiveness and anti-noise ability of the inversion method were verified by comparative experiments. The inversion results of fracture reservoir parameters show that the method has high applicability and accuracy in the inversion of porosity and fracture density; for the inversion results of fracture radius, the inversion accuracy of the method of the present invention in reservoirs containing small-scale fractures is significantly better than that in reservoirs containing large-scale fractures. This shows that the frequency-variant AVO inversion and parameter prediction method of orthogonal fracture reservoirs based on the Bayesian framework proposed in the present invention has good application prospects in fracture reservoir parameter estimation, especially when dealing with small-scale fractures, it shows higher accuracy and reliability.

[0173] Of course, the above description is not a limitation of the present invention, and the present invention is not limited to the above examples. Changes, modifications, additions or substitutions made by technicians in this technical field within the essential scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. A frequency-variable AVO inversion and parameter prediction method for orthogonal fracture reservoirs, characterized in that: The steps include: S1. Based on the Chapman model of multiple groups of fracture-pore media, the rock physics model of fluid-saturated horizontal-vertical orthogonal fracture media is derived and reconstructed to achieve rock physics modeling of fluid-saturated media containing two groups of orthogonal fractures; S2. Based on the anisotropic reflectivity algorithm and the rock physics model of fluid-saturated horizontal-vertical orthogonal fracture media, the frequency-variant AVO forward modeling of fractured reservoirs is realized, and the frequency-variant AVO reflection coefficient and synthetic seismic records of fractured reservoirs under the interface reflection model are calculated; S3. Driven by the forward modeling process of frequency-variable AVO, a frequency-variable AVO inversion theory based on the Bayesian framework is constructed. Using prestack seismic AVO data, the simulated annealing-particle swarm algorithm is introduced into the inversion process to realize frequency-variable AVO inversion based on the Bayesian inversion framework and predict reservoir parameters.

2. The method for frequency-variable AVO inversion and parameter prediction of orthogonal fracture reservoir according to claim 1, characterized in that: In step S1, the specific steps are as follows: Assume that there are two sets of orthogonal fractures f1 and f2 in the rock, and use the polarization angle θ m and azimuth Two angles characterize the distribution of cracks, m = 1, 2 represents two groups of cracks, [X1, X2] crack plane is the horizontal crack plane, [X1, X3] is the vertical crack plane, and it is defined that crack f1 is on the horizontal crack plane and crack f2 is on the vertical crack plane, that is, θ1 = 0°, θ2=90°, Then the normal direction of the crack surface is defined as: n 1 ＝(0,0,1)(1)； n 2 ＝(1,0,0)(2)； The frequency-dependent effective stiffness matrix of the porous equivalent medium with orthogonal cracks is derived as follows: In the formula, C ij (ω) is the equivalent elastic tensor of rock, is the elastic tensor of the isotropic background matrix, is the elastic correction term caused by porosity, is the elastic correction term caused by horizontal cracks, is the elastic correction term caused by vertical cracks; Elastic tensor of the isotropic background matrix Expressed by elastic parameters Lame coefficients λ and μ, the expression is: Elastic correction due to porosity The standard matrix form of is: In the formula, the expressions of d and e are: In the formula, is the porosity, v is the Poisson's ratio, is the pore stress tensor, p * is the fluid pressure in the pore; Pore ​​stress tensor and the fluid pressure in the pore, p * The expressions are: In the formula, is the porosity of group m fractures, τ m is the relaxation time scale parameter of m groups of cracks, is the critical pressure of m groups of cracks, K p is the bulk modulus of the pore, is the bulk modulus of m groups of cracks; Bulk modulus K of pores p , the bulk modulus of m groups of cracks and the critical pressure of group m fractures The expressions are: In the formula, κ f is the bulk modulus of the fluid; Elastic correction due to horizontal cracks and the elastic correction term due to vertical cracks The elastic parameter matrix of is characterized by the fluid pressure in the fracture and is expressed as: In the formula, the matrix elements The expression is: In the formula, r m is the aspect ratio of the crack, is the fluid pressure in group m fractures; m=1,θ1=0°, For the first group of fractures, the fluid pressure of the fractures arranged horizontally is expressed as: m=2,θ2=90°, For the second group of fractures, the fluid pressure of the fractures arranged vertically is expressed as: According to formula (3), the vertical P-wave velocity V of porous rock with orthogonal fractures saturated with fluid is obtained as follows: p (v), vertical shear wave velocity V s (ω) and the attenuation factor The expression is: In the formula, Re(C 33 (ω)) and Im(C 33 (ω)) represent the real part and imaginary part of the complex number respectively.

3. The frequency-variable AVO inversion and parameter prediction method for orthogonal fracture reservoirs according to claim 1, characterized in that: Step S2 specifically includes: The saturated fluid orthogonal fracture porous rock physics model is derived into a frequency-dependent anisotropic elastic model, and the equation is simplified by combining the anisotropic reflectivity method. After Fourier transformation, the expression of the anisotropic reflectivity method based on the differential equation system is obtained as follows: Where z is the depth; b is the displacement-stress vector; ω is the angular frequency; the matrix A is a 6×6 system matrix, which is composed of the elastic tensor C of the corresponding layer ij (ω), density ρ and horizontal slowness components in x and y directions; Based on the saturated fluid orthogonal fracture porous rock physics model, the upgoing wave reflection coefficient matrix R is re-derived using the eigenvalues ​​and eigenvectors of the system matrix A U , downlink wave reflection coefficient matrix R D , upgoing wave transmission coefficient matrix T U and the downlink transmission coefficient matrix T D , the expression is: Where, the matrix elements are the reflection responses in the frequency domain; When the energy is incident from top to bottom, the reflection and transmission coefficient matrix R of the anisotropic medium is expressed as: In the formula, R U , R D , T U , T D are the reflection coefficient matrix and transmission coefficient matrix of the upgoing and downgoing waves; The propagation matrix operator is introduced into the anisotropic reflectivity method and the reflection and transmission coefficient matrices are calculated through recursive equations. The reflection and transmission coefficient matrices of the i+1th layer are calculated by the i-th layer recursive relationship. The recursive relationship equation is: In the formula, r U 、r D ,t U ,t D The interface depth z=z i+1 The reflection matrix and transmission matrix at E U is the eigenvector of the upgoing wave reflection coefficient, E D is the eigenvector of the downlink wave reflection coefficient; E U 、E D The eigenvectors are: In the formula, h = z i+1 -z is the depth difference between the i-th layer and the i+1-th layer; λ is the eigenvalue of the corresponding system matrix A The frequency-varying AVO reflection coefficient matrix in the entire layered medium is obtained by iterative calculation using equations (37)-(40). The time domain prestack synthetic seismic record is obtained by inverse Fourier transform calculation of the frequency-varying AVO reflection coefficient and the seismic wavelet spectrum, and the formula is: In the formula, S PP is the pre-stack synthetic seismic record in the time domain; R PP (ω,θ) is the frequency-varying AVO reflection coefficient matrix, θ is the incident angle of the seismic wavelet, and F(ω) is the seismic wavelet spectrum.

4. The method for frequency-variable AVO inversion and parameter prediction of orthogonal fracture reservoir according to claim 1, characterized in that: In step S3, under the Bayesian inversion framework, the posterior probability formula of pre-stack seismic data for reservoir parameters is derived as follows: In the formula, m(ε f1 , ε f2 , Φ p , a f1 , a f2 ) is the model reservoir parameter; d is the model observation data, i.e., prestack seismic data; P(m) is the prior information of the model reservoir parameter; P(d|m)) is the likelihood function; ∫P(d|m))·P(m))dm is the probability normalization factor, which is a constant; The observed data d and the forward synthetic seismic record S are calculated by the calculation model PP The error between them, the objective function is constructed as: ΔE=∑||d-S PP ||2 (45); By exponentially transforming the error function, the likelihood function P(d|m(ε f1 , ε f2 , Φ p , a f1 , a f2 )): P(d|m(e f1 ,he f2 ,F p ,a f1 ,a f2 ))∝a·exp(-b,ΔE) (46); Where b is a constant; Obtaining prior information P(m(ε) from logging information f1 , ε f2 , Φ p , a f1 , a f2 ), using the model observation data as constraints, the reservoir parameter m(ε f1 , ε f2 , Φ p , a f1 , a f2 ), the objective function is constructed as: m(e f1 ,he f2 ,F p ,a f1 ,a f2 )=arg max[P(m(ε f1 ,he f2 ,F p ,a f1 ,a f2 )|d)] (47); The simulated annealing-particle swarm algorithm is used to solve the objective function formula (47). By searching for the parameter m(ε f1 , ε f2 , Φ p , a f1 , a f2 ) is the optimal solution, so that the posterior probability P(m(ε f1 , ε f2 , Φ p , a f1 , a f2 )|d)maximum; end the calculation and output the inversion results.

Images