A method for evaluating fracture development in a tight sand gas reservoir on a well logging scale

Abstract

The application discloses a kind of compact sandstone gas reservoir logging scale fracture development evaluation method, it is related to the field of earth exploration physics, including: based on complex pore structure high-angle fracture type rock physics model and actual logging data, two groups of fracture density parameters in well are obtained by using simulated annealing particle swarm algorithm inversion;Respectively calculate the fracture density identification factor under different rotation angles;Based on the correlation coefficient between Poisson impedance and fracture density under different rotation angles, determine the optimal rotation angle when the correlation coefficient reaches maximum;Combining longitudinal wave impedance and shear wave impedance data, calculate the fracture density identification factor under the optimal rotation angle;Fracture density identification factor under the optimal rotation angle is applied to actual logging data, and fracture development is characterized.The method can effectively distinguish and quantify the multiple sets of high-angle fractures developed in compact sandstone reservoir, significantly improve the prediction accuracy and reliability of fracture reservoir, and provide key basis for reservoir sweet spot prediction and well deployment.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of earth exploration physics technology, and particularly relates to a method for evaluating fracture development in a tight sandstone gas reservoir on a logging scale. BACKGROUND

[0002] As an important part of unconventional oil and gas resources, the position of tight sandstone gas in the global energy structure is increasingly prominent. Compared with conventional reservoirs, tight sandstone reservoirs usually have the inherent characteristics of low porosity and low permeability, and their economic exploitation value largely depends on the natural fracture system developed in the reservoir. These fractures are not only important spaces for oil and gas accumulation, but also key channels for fluid percolation.

[0003] Rock physics model is a bridge connecting the microstructure of reservoir (mineral, pore, fracture, fluid) and the macroscopic elastic response (velocity, density, impedance). Traditional single rock physics model usually assumes that the fracture is dry or single group, and ignores the complex pore structure of the matrix. Tight sandstone reservoirs are essentially typical multi-scale pore-fracture dual media, containing matrix pores from nanoscale to microscale, millimeter-scale microfractures, and even centimeter-scale large fractures. In addition, the reservoir is usually in a partially saturated state, i.e., natural gas and formation water exist in the pore space at the same time. Under such complex physical conditions, any single theoretical model is difficult to accurately describe the elastic properties. In terms of parameter inversion, estimating multiple fracture parameters from limited observation data is a typical underdetermined and nonlinear optimization problem. The traditional linear inversion method is heavily dependent on the initial model and is prone to fall into a local optimal solution, and cannot obtain a stable and reliable global optimal solution. In terms of fracture identification factors, many empirical or semi-empirical indicators based on conventional logging curves have been proposed, such as resistivity-acoustic travel time crossplot method, P-S wave velocity ratio, etc. These methods can indicate fracture development intervals to some extent, but they generally have multiple solutions, cannot effectively distinguish multiple groups of fractures, and are sensitive to fluid response, with similar characteristics in gas layers and water layers, and unstable application effect.

[0004] Therefore, there is an urgent need in the art to develop a new fracture evaluation method for tight sandstone gas reservoirs that can overcome the above limitations. SUMMARY

[0005] To address the aforementioned issues, this invention discloses a logging-scale fracture development evaluation method for tight sandstone gas reservoirs. First, this method constructs a multi-scale partially saturated petrophysical model that comprehensively considers mineral matrices such as quartz and clay, matrix porosity, intergranular porosity, microfractures, and two sets of high-angle orthogonal vertical fractures. This model accurately describes the complex pore structure and fluid partial saturation effect of tight sandstone, achieving precise equivalence from the mineral matrix to the complex fracture-pore medium. Second, using the constructed model as a forward modeling engine, and employing a simulated annealing particle swarm optimization algorithm, with measured logging P-wave velocities and S-wave velocities as constraints, the density parameters of the two sets of fractures are inverted, achieving quantitative decoupling of the complex fracture system. Finally, the method innovatively introduces and modifies Poisson impedance theory to construct a fracture density identification factor that is sensitive to fracture density but relatively insensitive to fluid effects. This factor effectively eliminates the interference of fluid effects and, through cross-correlation matching of elastic impedance data at different angles, greatly suppresses the contribution of pore fluid variations, thereby highlighting the sensitive response to fracture density, especially the difference in density between the two sets of fractures. The method of this invention can effectively distinguish and quantify multiple sets of high-angle fractures developed in tight sandstone reservoirs, significantly improving the accuracy and reliability of fractured reservoir prediction, and providing key basis for reservoir sweet spot prediction and well location deployment.

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

[0007] A method for evaluating fracture development at the logging scale in tight sandstone gas reservoirs includes the following steps:

[0008] s1. Based on a rock physics model of high-angle fractures with complex pore structures and actual well logging data, two sets of fracture density parameters were obtained by inverting the model using a simulated annealing particle swarm optimization algorithm, with actual well logging P-wave velocity and S-wave velocity as constraints. , ;

[0009] s2. Based on well-logging-scale P-wave impedance and S-wave impedance data, this paper innovatively introduces and improves Poisson impedance theory, proposing a fracture density identification factor; based on well-logging-scale P-wave impedance... and transverse wave impedance Data, calculate different rotation angles respectively Lower Crack Density Identification Factor and ;

[0010] s3. Based on formula Calculate Poisson impedance and crack density at different rotation angles , The correlation coefficient between them is used to determine the optimal rotation angle when the correlation coefficient reaches its maximum value. and ;

[0011] s4. Calculate the fracture density identification factor under the optimal rotation angle combined with the logging scale P-wave impedance and S-wave impedance data and ;

[0012] s5. Apply the fracture density identification factor under the optimal rotation angle to the actual logging data to characterize the fracture development.

[0013] Optionally, in step s1, the construction of the complex pore structure high-angle fracture type rock physics model comprises:

[0014] s11. Calculate the equivalent modulus of the solid matrix

[0015] The equivalent bulk modulus of the solid matrix is calculated using the VRH average theory and the equivalent shear modulus of the solid matrix , whose expressions are as follows:

[0016] ;

[0017] ;

[0018] wherein, is the inclusion phase, is the number of inclusion phase types, is the volume fraction of the inclusion phase, is the bulk modulus of the inclusion, is the shear modulus of the inclusion, is the equivalent bulk modulus of the solid matrix, is the equivalent shear modulus of the solid matrix;

[0019] s12. Construct a matrix dry skeleton containing multiple pores

[0020] The micro-pores are included in the solid matrix using the K-T model. For a solid matrix dry skeleton containing randomly distributed micro-pores, the equivalent bulk modulus and the shear modulus of the solid matrix dry skeleton are expressed as follows:

[0021] ;

[0022] ;

[0023] ;

[0024] ;

[0025] wherein, is the porosity of the micro-pores, and geometric factor of inter-particle pores, equivalent bulk modulus of solid matrix dry skeleton containing randomly distributed inter-particle pores, equivalent shear modulus of solid matrix dry skeleton containing randomly distributed inter-particle pores;

[0026] The differential effective medium theory is used to gradually add inter-particle pores of different shapes to the solid matrix dry skeleton containing randomly distributed inter-particle pores to construct a matrix dry skeleton containing multiple pores, and the equivalent bulk modulus and the shear modulus of the matrix dry skeleton containing multiple pores are expressed as follows:

[0027] ;

[0028] ;

[0029] wherein, porosity of inter-particle pores, and geometric factor of inter-particle pores, equivalent bulk modulus of matrix dry skeleton containing multiple pores, equivalent shear modulus of matrix dry skeleton containing multiple pores;

[0030] s13. Constructing a dry skeleton with developed micro-cracks

[0031] The Hudson model is used to add directionally arranged micro-cracks to the matrix dry skeleton containing multiple pores to obtain a dry skeleton with developed micro-cracks having transverse isotropy, and the equivalent stiffness matrix is as follows:

[0032] ;

[0033] wherein, isotropic stiffness matrix of matrix dry skeleton containing multiple pores, first-order micro-crack perturbation term, second-order micro-crack perturbation term;

[0034] For the convenience of model calculation, only the first-order micro-crack perturbation term given by the Hudson model is calculated, and the specific expression is as follows:

[0035] ;

[0036] ;

[0037] ;

[0038] ;

[0039] ;

[0040] where, and is the Lame constant of the porous matrix dry skeleton, is the crack density of the micro-cracks, and is a parameter depending on the crack aspect ratio and the background medium modulus;

[0041] s14. Constructing a dry rock skeleton with cracks

[0042] Using the Schoenberg linear slip theory, on the basis of the dry skeleton with micro-cracks, two groups of macroscopic high-angle vertical cracks are introduced, and the two groups of cracks are orthogonal in azimuth. A dry rock skeleton with cracks is constructed. The equivalent stiffness matrix of the dry rock skeleton with cracks is obtained by the stiffness matrix of the background medium and the crack compliance matrix. The model assumes that the normal direction of the first group of cracks is along the x-axis, and the normal direction of the second group of cracks is along the y-axis. The normal and tangential compliance matrix expressions are as follows:

[0043] ;

[0044] ;

[0045] ;

[0046] ;

[0047] where, is the normal compliance of the first group of cracks, is the tangential compliance of the first group of cracks, is the normal compliance of the second group of cracks, is the tangential compliance of the second group of cracks, is the crack density of the first group of cracks, is the crack density of the second group of cracks, is the crack aspect ratio of the first group of cracks, is the crack aspect ratio of the second group of cracks, is the Poisson's ratio function;

[0048] The equivalent stiffness matrix expression of the dry rock skeleton with cracks is as follows:

[0049] ;

[0050] where, is the crack compliance matrix, which is a specific form of matrix composed of , , , ;

[0051] s15. Constructing a complex-porosity high-angle fracture rock physics model

[0052] By introducing the microscale and mesoscale dual fluid flow mechanisms to simulate the dispersion and attenuation effects under partial saturation conditions, filling the mixture of partially saturated fluids into the dry rock skeleton containing fractures, a partially saturated tight sandstone rock, i.e., a complex-porosity high-angle fracture rock physics model, is constructed, and the equivalent modulus expression of the model is as follows:

[0053] ;

[0054] ;

[0055] ;

[0056] wherein, is the frequency-dependent bulk modulus of the partially saturated tight sandstone rock, is the bulk modulus of the dry rock skeleton containing fractures, is the frequency-dependent shear modulus of the partially saturated tight sandstone rock, is the rock porosity, is the microfracture porosity, is the first group of fracture porosities, is the second group of fracture porosities, is the angular frequency, is the bulk modulus of the mixed fluid, is the imaginary unit, is the characteristic frequency determined by the fluid diffusion coefficient;

[0057] The theoretical velocity of the partially saturated tight sandstone rock, i.e., the longitudinal wave velocity and the transverse wave velocity of the model, is calculated as follows:

[0058] ;

[0059] ;

[0060] ;

[0061] wherein, is the longitudinal wave velocity of the model, is the transverse wave velocity of the model, is the solid matrix density, is the density of water, is the density of gas, is the water saturation, is the gas saturation.

[0062] Optionally, in step s1, the two groups of fracture density parameters in the well are obtained by inversion 、 The step of the step s1, specifically comprises: matching the actual shear wave velocity curve and the longitudinal wave velocity curve in the logging data with the theoretical shear wave velocity curve and the longitudinal wave velocity curve calculated by the model, constructing the objective function of inversion, and defining the objective function As the error norm between the model prediction value and the measured value, and the expression is as follows:

[0063] ;

[0064] Wherein, The model parameter vector to be inverted is The measured longitudinal wave velocity in the well is The measured shear wave velocity in the well is And The weight coefficient is The regularization term is used to introduce prior information.

[0065] Optionally, in step s2, the fracture density identification factor is calculated based on the Poisson impedance method, and the expression is as follows:

[0066] ;

[0067] ;

[0068] Wherein, The rotation angle is determined by optimizing the petrophysical analysis, The longitudinal wave impedance is The shear wave impedance is.

[0069] The beneficial effects of the present application are as follows:

[0070] (1) Comprehensive consideration of reservoir complexity. The petrophysical model established by the present application fully considers the characteristics of multiple mineral components, multiple scale pore structures, two groups of high-angle orthogonal vertical fractures and partially saturated fluids in the dense sandstone gas reservoir, and can better reflect the complex physical characteristics of the actual reservoir than the traditional model, thereby laying a solid foundation for improving the fracture prediction accuracy.

[0071] (2) Strong decoupling ability of inversion method. The simulated annealing particle swarm optimization algorithm adopted has global search ability and local fine optimization ability, effectively overcomes the local optimal problem in nonlinear inversion, ensures the stability and reliability of the inversion result under the complex model, successfully realizes the quantitative and synchronous inversion of the two groups of fracture densities, and solves the problem of distinguishing multiple groups of fractures.

[0072] (3) Effectively decoupling multiple sets of fractures and fluid effects. One of the core innovations of the present application is to propose a fracture density identification factor, which, through mathematical transformation, cleverly amplifies the response of the fractures while suppressing the response of the fluid. Through inversion, the fracture density of the two sets of fractures is directly obtained, realizing the separate quantification of the two sets of fractures, breaking through the limitation of traditional methods that can only identify the dominant fracture set; cleverly enhancing the sensitivity to fractures while suppressing the sensitivity to fluid.

[0073] (4) The evaluation results are intuitive and quantitative. The final output results not only include two fracture density curves, but also provide a fracture density identification factor curve. The fracture density identification factor is highly correlated with the fracture density, making the identification of fracture development degree very intuitive and quantitative, greatly facilitating the use of geology interpreters and being conducive to quickly locking the "sweet spot" interval. BRIEF DESCRIPTION OF DRAWINGS

[0074] Figure 1 A flow chart of the present application for evaluating fracture development in a tight sand gas reservoir on a logging scale.

[0075] Figure 2 A comprehensive logging curve of the Y well in the Shihexizi group tight sand gas reservoir shown in the embodiment of the present application.

[0076] Figure 3 Fracture density and simulated velocity curve results based on a rock physics model inversion of the Y well in the study area shown in the embodiment of the present application.

[0077] Figure 4 Correlation between the fracture density identification factor and the elastic impedance of the Y well in the study area shown in the embodiment of the present application, wherein (a) is the correlation coefficient changes with the rotation angle, and (b) is the correlation coefficient changes with the rotation angle.

[0078] Figure 5 The intersection of and with the color scale of and the intersection of and the fracture density identification factor at the optimal rotation angle, wherein (a) is the intersection of the compressional wave impedance and the shear wave impedance (represented by the first set of fracture density), and (b) is the intersection of the compressional wave impedance and the first set of fracture density identification factor .

[0079] Figure 6 The intersection of and with the color scale of and the intersection of The intersection of the fracture density identification factor at the optimal rotation angle with the intersection of the fracture density identification factor at the optimal rotation angle , wherein (a) is the intersection of the P-wave impedance and S-wave impedance (expressed in the second set of fracture density), and (b) is the intersection of the P-wave impedance and the second set of fracture density identification factor. .

[0080] Figure 7 The fracture density identification factor curve of the research area Y well at the optimal rotation angle is shown in the embodiment of the present application. DETAILED DESCRIPTION

[0081] In order to make the purpose, technical scheme and advantages of the embodiments of the present application clearer, the technical scheme in the embodiments of the present application will be described clearly and completely below in conjunction with the drawings in the embodiments of the present application. Obviously, the described embodiments are part of the embodiments of the present application, rather than all the embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by a person of ordinary skill in the art without creative labor fall within the protection scope of the present application. Therefore, the following detailed description of the embodiments of the present application provided in the drawings is not intended to limit the scope of the claimed present application, but only represents selected embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by a person of ordinary skill in the art without creative labor fall within the protection scope of the present application.

[0082] In view of the deficiencies of the prior art, the present application proposes a new tight sandstone gas reservoir logging scale fracture development evaluation method, and the specific purposes include: (1) building a rock physics model with complete physical mechanism, which can finely represent the multi-scale pore structure, two sets of high-angle orthogonal vertical fractures and partial saturated fluid characteristics of the tight sandstone reservoir, laying a theoretical foundation for accurately predicting rock elastic properties; (2) providing a stable and efficient fracture parameter inversion method, which solves the problem that the traditional linear inversion method is easy to fall into local optimal solution under a complex nonlinear rock physics model by introducing a simulated annealing particle swarm hybrid optimization algorithm, and realizes global optimal inversion of two sets of fracture density parameters; (3) inventing a new fracture density identification factor, which is not sensitive to fluid changes and is highly sensitive to fracture density, especially the difference between the two sets of fracture density, so as to realize direct, continuous and quantitative evaluation of fracture development degree by using only conventional logging data, and provide a new and more reliable quantitative index for reservoir evaluation; (4) ultimately forming a complete technical process from physical modeling, parameter inversion to comprehensive evaluation, and providing reliable geophysical basis for "sweet spot" prediction, well deployment and productivity evaluation of tight sandstone gas reservoirs.

[0083] A tight sandstone gas reservoir logging scale fracture development evaluation method, comprising the following steps:

[0084] s1. Based on the complex pore structure high-angle fracture type rock physics model and the actual logging data, the simulated annealing particle swarm optimization algorithm is used to obtain two groups of fracture density parameters by inversion taking the actual logging longitudinal wave velocity and transverse wave velocity as constraints 、 .

[0085] 1. Complex pore structure high-angle fracture type rock physics modeling

[0086] As an important unconventional natural gas resource, tight sandstone gas occupies an increasingly important position in the global energy structure. By comprehensively utilizing core, thin section, scanning electron microscope and imaging logging data, the tight sandstone reservoir is finely analyzed. The main minerals of the reservoir matrix are quartz, feldspar and clay, etc. The reservoir usually has inherent characteristics of low porosity and low permeability, and has typical multi-scale pore structure, including matrix pore, intergranular micro-pore, diagenetic micro-fracture and macro natural tectonic fracture. The development degree of the natural fracture system is a key geological factor controlling high yield and stable yield of the reservoir, especially the high-angle fractures. They are not only the dominant channel for oil and gas migration, but also the decisive factor for forming effective percolation network and greatly improving the permeability of the reservoir.

[0087] The present application aims to establish a rock physics model which can fully reflect the complex characteristics of the tight sandstone gas reservoir. The model comprehensively considers the multi-mineral components of the reservoir rock, multi-scale pore structure (including matrix pore, micro-fracture and macro fracture), two groups of high-angle orthogonal vertical fractures and the distribution characteristics of partially saturated fluid. By introducing the complex multi-scale pore system and the double fluid mechanism, the model can more accurately describe the complex physical characteristics of the actual tight sandstone gas reservoir, and lay a model foundation for accurately inverting the fracture parameters. The construction of the complex pore structure high-angle fracture type rock physics model includes:

[0088] s11. Calculating the equivalent modulus of the solid matrix

[0089] The tight sandstone is usually mainly composed of quartz and clay, and may contain a small amount of other minerals. The Voigt-Reuss-Hill (VRH) average theory is used to calculate the equivalent bulk modulus of the solid matrix and the equivalent shear modulus of the solid matrix , and the expressions are as follows:

[0090] (1);

[0091] (2);

[0092] wherein, is the included phase, is the phase class number of the included phase, is the volume fraction of the included phase, the bulk modulus of the inclusion, the shear modulus of the inclusion, the equivalent bulk modulus of the solid matrix, the equivalent shear modulus of the solid matrix.

[0093] s12. Constructing a matrix dry skeleton containing multiple pores

[0094] The pore system of tight sandstones is complex, including micropores associated with clay minerals and relatively larger intergranular pores. The influence of the ubiquitous micropores in the matrix on the elastic properties is similar to that of spherical pores. The Kuster-Toksöz (K-T) model is used to include the micropores into the solid matrix. For a solid matrix dry skeleton containing randomly distributed micropores, the equivalent bulk modulus and the shear modulus are expressed as follows:

[0095] (3);

[0096] (4);

[0097] (5);

[0098] (6);

[0099] where, is the porosity of the micropores, and is the geometric factor of the micropores, is the equivalent bulk modulus of the solid matrix dry skeleton containing randomly distributed micropores, is the equivalent shear modulus of the solid matrix dry skeleton containing randomly distributed micropores;

[0100] Tight sandstones often contain dispersed relatively rigid intergranular pores, which are usually non-spherical in shape. The differential effective medium theory (DEM) can flexibly handle pores of different aspect ratios. The differential effective medium theory (DEM) is used to gradually add intergranular pores of different shapes to the solid matrix dry skeleton containing randomly distributed micropores to construct a matrix dry skeleton containing multiple pores. The equivalent bulk modulus and the shear modulus are expressed as follows:

[0101] (7);

[0102] (8);

[0103] where, is the porosity of the intergranular pores, and is the geometric factor of intergranular pores, is the dry skeleton equivalent bulk modulus of the matrix containing multiple pores, is the dry skeleton equivalent shear modulus of the matrix containing multiple pores; the calculation process is from the pure matrix , , Start with the solid matrix dry skeleton containing randomly distributed spherical pores, gradually add intergranular pores with different aspect ratios to the background, until the porosity reaches , and finally get and , which are the equivalent modulus of the matrix dry skeleton containing multiple pores.

[0104] s13. Constructing a dry skeleton with microfracture development

[0105] Microfractures usually develop in tight sandstone reservoirs, with small size, good connectivity, usually considered as flat ellipsoidal inclusions with small aspect ratio, and have certain dominant orientation arrangement direction. The Hudson model is used to add directional microfractures to the matrix dry skeleton containing multiple pores to obtain a dry skeleton with transversely isotropic microfracture development, and its equivalent stiffness matrix is as follows:

[0106] (9);

[0107] wherein, is the isotropic stiffness matrix of the matrix dry skeleton containing multiple pores, is the first-order microfracture perturbation term, is the second-order microfracture perturbation term;

[0108] For convenience of model calculation, only the first-order microfracture perturbation term given by the Hudson model is calculated, and its specific expression is as follows:

[0109] (10);

[0110] (11);

[0111] (12);

[0112] (13);

[0113] (14);

[0114] wherein, and are the Lame constants of the matrix dry skeleton containing multiple pores, The crack density is the density of the microcracks. and It depends on the crack aspect ratio and the modulus of the background medium.

[0115] s14. Construct a dry rock skeleton containing cracks

[0116] For the two sets of high-angle orthogonal vertical fracture systems commonly found in tight sandstone gas reservoirs, orthogonal anisotropic media is the geological model that best matches the actual strata. Schoenberg linear slip theory treats fractures as infinitely thin planes with linear slip boundary conditions, representing the influence of fractures as an additional flexibility of the stiffness matrix. For the orthogonal vertical fracture system, the total fracture flexibility matrix is ​​the sum of the two sets of fracture flexibility matrices. Using Schoenberg linear slip theory, based on a dry framework with microfractures, two sets of macroscopic high-angle vertical fractures are introduced, and the two sets of fractures are orthogonal to each other in orientation, constructing a fractured dry rock framework. The equivalent stiffness matrix of the fractured dry rock framework is obtained from the stiffness matrix of the background medium and the fracture flexibility matrix. The model assumes that the normal direction of the first set of fractures is along the x-axis, and the normal direction of the second set of fractures is along the y-axis. The expressions for their normal and tangential flexibility matrices are as follows:

[0117] (15);

[0118] (16);

[0119] (17);

[0120] (18);

[0121] in, For the first group of crack normal compliance, For the first set of crack tangential compliance, For the second group of crack normal compliance, For the second group of crack tangential compliance. The crack density of the first group of cracks. The crack density of the second group of cracks. The aspect ratio of the first group of cracks. The aspect ratio of the second group of cracks. It is the Poisson ratio function;

[0122] The equivalent stiffness matrix expression for a cracked, dry rock skeleton is as follows:

[0123] (19);

[0124] in, The crack compliance matrix is ​​derived from... 、 、 、 constituted in a particular form of matrix.

[0125] s15. Constructing a complex-porosity high-angle fracture rock physics model

[0126] Tight sand gas reservoirs are usually in a partially saturated state, i.e. containing gas and liquid in the pore space simultaneously, and the velocity dispersion and attenuation caused by local pressure imbalance due to uneven fluid distribution. Usually for microscale fluid flow, an equivalent tubular pore model is used, and for mesoscale fluid flow, a White spherical patch saturation model is used. By introducing microscale and mesoscale double fluid flow mechanisms to simulate the dispersion and attenuation effects under partial saturation conditions, filling the partially saturated fluid mixture into the dry rock skeleton containing fractures, a partially saturated fluid tight sandstone rock, i.e. a complex-porosity high-angle fracture rock physics model, is constructed, and the equivalent modulus expression of the model is as follows:

[0127] (20);

[0128] (21);

[0129] (22);

[0130] wherein, is the frequency-dependent bulk modulus of the partially saturated tight sandstone rock, is the bulk modulus of the dry rock skeleton containing fractures, is the frequency-dependent shear modulus of the partially saturated tight sandstone rock, is the rock porosity, is the microfracture porosity, is the first group of fracture porosity, is the second group of fracture porosity, is the angular frequency, is the bulk modulus of the mixed fluid, is the imaginary unit, is the characteristic frequency determined by the fluid diffusion coefficient;

[0131] The theoretical velocity of the partially saturated tight sandstone rock is thus calculated, i.e. the longitudinal wave velocity and the transverse wave velocity of the model, and the expressions are as follows:

[0132] (23);

[0133] (24);

[0134] (25);

[0135] where, is the model P-wave velocity, is the model S-wave velocity, is the solid matrix density, is the water density, is the gas density, is the water saturation, is the gas saturation.

[0136] 2. Fracture density inversion based on simulated annealing particle swarm optimization algorithm

[0137] In geophysical exploration, inversion of reservoir parameters based on rock physics model and logging data is the core step of accurate characterization of reservoir. However, the inversion process has high nonlinearity, multi-parameter coupling and multiple extremum, and the traditional optimization algorithm is difficult to effectively solve. Simulated annealing particle swarm optimization algorithm (SA-PSO) as an intelligent hybrid optimization algorithm shows significant advantages and effects in this process.

[0138] The fracture system of tight sandstone gas reservoir is complex, and the fracture parameters have obvious multi-scale and multi-directional distribution characteristics. In order to realize the quantitative identification of multiple sets of fractures, the present application proposes a fracture density inversion method based on simulated annealing particle swarm optimization algorithm. The core is to use the measured logging data as the constraint based on the complex pore structure high-angle fracture type rock physics model, use the constructed model as the forward engine, use the simulated annealing particle swarm optimization algorithm, jointly invert the P-wave and S-wave velocities and the measured logging data, and realize the accurate prediction of the densities of two sets of fractures. The actual S-wave velocity curve and P-wave velocity curve in the logging data are matched with the theoretical S-wave velocity curve and P-wave velocity curve calculated by the model, the objective function of inversion is constructed, and the objective function is defined as is the error norm between the model predicted value and the measured value, and the expression is as follows:

[0139] (26);

[0140] where, is the model parameter vector to be inverted, is the measured P-wave velocity in the well, is the measured S-wave velocity in the well, and are weight coefficients, is a regularization term, which is used to introduce prior information such as non-negative fracture density, spatial smoothing constraint, etc., to improve the stability of inversion.

[0141] Since the objective function is relative to the model parameter The nonlinear, multimodal, the standard particle swarm algorithm is easy to fall into local extremum, and convergence precision is insufficient for multi-parameter, nonlinear strong optimization problem. The simulated annealing algorithm is introduced, the standard particle swarm algorithm is improved, and the hybrid optimization algorithm is formed, the hybrid optimization algorithm combines the fast convergence of the particle swarm algorithm and the probability jump characteristics of the simulated annealing algorithm, and the global optimal solution can be effectively searched. The specific steps of the algorithm are as follows:

[0142] (1) initialize particle swarm and simulated annealing parameters

[0143] N particles are randomly generated in the feasible region of the parameter space, to form an initial population, the position and speed of each particle are randomly initialized, the position of each particle represents the candidate solution of two groups of crack density, meanwhile, the maximum iteration number is set, and the initial temperature and temperature reduction coefficient of simulated annealing are initialized.

[0144] (2) particle fitness value initialization calculation and evaluation

[0145] For each particle, according to the crack density combination represented by it, the established rock physical model is used to simulate the longitudinal wave velocity and transverse wave velocity, and is matched with the actual logging data, the objective function value is calculated, the fitness value of each particle is evaluated, the individual optimal position of each particle is initialized, and the global optimal position of the whole population is found.

[0146] (3) particle velocity and position update

[0147] Iterative calculation is carried out, and the velocity and position of each particle are updated according to the individual and group optimal position.

[0148] (4) individual and global optimal update

[0149] For each particle, the current fitness of each particle is compared with the historical optimal fitness, the particle individual optimal solution and the corresponding particle position are updated, and the particle with the best fitness in the current population is found as the global optimal solution.

[0150] (5) simulated annealing disturbance

[0151] The global optimal solution is disturbed by simulated annealing, a new global optimal solution is generated, and the fitness difference is calculated. If the fitness difference is less than 0, the new global optimal solution is accepted; if the fitness difference is greater than or equal to 0, the new global optimal solution is accepted with a temperature reduction probability.

[0152] (6) temperature reduction and iteration

[0153] According to the temperature reduction plan, the temperature is reduced, and whether the termination condition (reaching the maximum iteration number or the fitness being small enough) is met is checked, if not, the particle velocity and position updating operation is carried out again.

[0154] (7) Convergence judgment

[0155] When the maximum number of iterations is reached or the objective function value converges to a predetermined threshold, the iteration is stopped; the final global optimal solution is output, that is, the two sets of fracture densities obtained by inversion. By performing the above inversion point by point (each depth point), finally two fracture density curves varying with depth are obtained.

[0156] s2. Based on the logging scale P-wave impedance and S-wave impedance data, the Poisson impedance theory is innovatively introduced and improved, and a fracture density identification factor is proposed. According to the logging scale P-wave impedance and S-wave impedance data, the fracture density identification factors and under different rotation angles are calculated respectively.

[0157] Although the fracture density is obtained by inversion based on the rock physical model and logging data, there are obvious limitations under the coupling condition of complex fracture system and fluid, mainly manifested as insufficient resolution ability for multiple groups of fractures, and mutual interference between fracture parameters and fluid effects, which reduces the contribution differentiation degree of elastic parameters to both, thereby limiting the fine characterization of the reservoir fracture system. It is urgent to have a more intuitive and more stable attribute to directly indicate the fracture development strength, and this attribute is preferably used for the well section without complex inversion or as the attribute target of seismic inversion. The present application further constructs a new type of fracture density identification factor based on the Poisson impedance theory. The factor is sensitive to the change of fracture density, relatively insensitive to the change of fluid and lithology, further highlights the fracture information, suppresses the residual fluid and lithology interference, and provides a more intuitive and reliable fracture development degree evaluation method.

[0158] P-wave impedance and S-wave impedance are important parameters for characterizing rock elastic properties, and it is difficult to accurately distinguish fracture development characteristics when they are used alone. Based on the proposed Poisson impedance theory, the present application constructs an identification factor for identifying the density of two groups of high-angle orthogonal vertical fractures in combination with the rock physical properties of the tight sandstone reservoir. The Poisson impedance (PI) method realizes coordinate rotation transformation through linear combination of P-wave and S-wave impedance, thereby enhancing the sensitivity to the change of specific rock physical properties. In view of the lack of azimuth anisotropy azimuth seismic data in the actual work area, the following fracture density identification factors are proposed, which can effectively characterize the development of high-angle orthogonal vertical fractures. The fracture density identification factor is calculated based on the Poisson impedance method, and the expression is as follows:

[0159] (27);

[0160] (28);

[0161] wherein, is the rotation angle determined by the optimization of petrophysical analysis, is the P-wave impedance, is the S-wave impedance.

[0162] s3. Calculate the correlation coefficient between the Poisson impedance and the fracture density under different rotation angles based on the formula , , , , .

[0163] s4. Calculate the fracture density identification factor under the optimal rotation angle in combination with the logging scale P-wave impedance and S-wave impedance data , .

[0164] s5. Apply the fracture density identification factor under the optimal rotation angle to the actual logging data to characterize the fracture development.

[0165] Application example

[0166] The application example of the present application carries out research and experiment on the logging scale fracture development of a tight sand gas reservoir in a region with a northeast-southwest overall trend. As a typical tectonic-sedimentary transition zone, the region has important research significance in terms of fracture development, sedimentary response and oil and gas accumulation. The Y well of the Shihezi Formation tight sand gas in the research region is selected as the research well site, and the reservoir of the well site has the characteristics of “low porosity, low permeability and strong heterogeneity”, which is a typical tight gas reservoir. Tectonic fractures and tectonic controlled natural fracture networks are widely developed in the formation, forming effective gas migration channels and being the key control factors for the formation of complex fracture networks in the process of hydraulic fracturing. Figure 2 is the logging curve of the Y well, including gamma ray GR, P-wave velocity , S-wave velocity , density , porosity and water saturation . The logging data shows the characteristics of “three lows and one high”: low gamma, low density, low water saturation and high velocity, indicating that the reservoir in the research region has good geological response characteristics of tight sand gas reservoir.

[0167] The fracture density inversion method based on the simulated annealing particle swarm algorithm proposed in the present application is used to invert the logging data of the Y well point by point (each depth point), and two fracture density curves varying with depth, model simulated S-wave velocity curve and P-wave velocity curve are obtained, as shown in Figure 3As shown, the longitudinal wave velocity obtained by the target interval inversion has a high degree of fitting with the measured longitudinal wave velocity in the well, the simulated transverse wave velocity has a good fitting with the measured transverse wave velocity in the well, and the fracture density curve presents obvious zonation characteristics at different depth intervals and is consistent with the changes of lithology, porosity and fluid-containing properties. These results not only verify the accuracy and applicability of the method in the quantitative prediction of fracture density in the tight sandstone reservoir, but also provide reliable technical support for the parameter characterization and development and utilization of the fractured reservoir.

[0168] By calculating the correlation coefficient of Poisson impedance and fracture density, the optimal rotation angles are determined as and , as shown in Figure 4 (a) and (b). Based on the optimal rotation angles, the fracture density identification factor and the Poisson impedance cross-plot capable of effectively characterizing the fracture density are calculated, as shown in Figure 5 (a), (b) and Figure 6 (a), (b). and The cross-plot of and the fracture density identification factor has weak differentiation of different fracture densities and no ability to identify the fracture density alone, while in the cross-plot of and the fracture density identification factor , different fracture density regions show obvious separation, and high fracture density (orange) corresponds to the low value region, further verifying the superiority of the parameter in identifying the fracture development degree. The results show that the fracture density identification factor parameter can effectively characterize the fracture development degree, and

[0169] can be used as a parameter to characterize the fracture development characteristics of the tight sandstone gas reservoir. Figure 7 The fracture density identification factor curve under the optimal rotation angle of Y well. The results show that in the target interval, the fracture density identification factor has a significant negative correlation with the fracture density, that is, the high fracture density region corresponds to the low fracture density identification factor value, and this corresponding relationship verifies the effectiveness of the fracture density identification factor in characterizing the fracture development of the tight reservoir.

[0170] Of course, the above description is not a limitation on the present application, and the present application is not limited to the above examples. Changes, modifications, additions or substitutions made by those skilled in the art within the essential scope of the present application should also be within the protection scope of the present application.

Claims

1. A method for evaluating fracture development at a well logging scale in a tight sand gas reservoir, characterized in that, The method comprises the following steps: s1. Based on the complex pore structure high angle fracture type rock physics model and the actual logging data, the simulated annealing particle swarm optimization algorithm is used to obtain two groups of fracture density parameters by inverting the actual logging longitudinal wave velocity and shear wave velocity as constraints 、 ; s2. Based on Poisson impedance theory, according to the longitudinal wave impedance at the logging scale and transverse wave impedance Data, calculate different rotation angles respectively Lower Crack Density Identification Factor and ; s3. Calculate the correlation coefficient between the Poisson impedance and the fracture density under different rotation angles, and determine the optimal rotation angle when the correlation coefficient reaches the maximum value 、 and ;​ s4. In combination with the well logging scale P-wave impedance and S-wave impedance data, the fracture density identification factor under the optimal rotation angle is calculated and ; s5. The fracture density identification factor at the optimal rotation angle is applied to the actual logging data to characterize the fracture development; In step S1, the construction of the complex pore structure high-angle fracture type rock physical model comprises: The equivalent modulus of the solid matrix is calculated by using the VRH average theory; The micro-pores are included in the solid matrix by using the K-T model to obtain a solid matrix dry skeleton containing randomly distributed micro-pores, and the differential equivalent medium theory is used to gradually add intergranular pores of different shapes to the solid matrix dry skeleton containing randomly distributed micro-pores to construct a matrix dry skeleton containing multiple pores; The micro-fractures in the directional arrangement are added to the matrix dry skeleton containing multiple pores to obtain a dry skeleton with micro-fracture development in the transverse isotropy; On the basis of the dry skeleton with micro-fracture development, two groups of macro high-angle vertical fractures are introduced to construct a dry rock skeleton containing fractures; The partially saturated fluid mixture is filled into the dry rock skeleton containing fractures to obtain a complex pore structure high-angle fracture type rock physical model; In step S2, the fracture density identification factor is calculated based on the Poisson impedance method, and the expression is as follows: ； ； wherein, is the rotation angle determined by petrophysical analysis optimization, is the P-wave impedance, is the S-wave impedance.

2. The method of claim 1, wherein the method is used to evaluate the development of fractures in a tight sand gas reservoir on a well logging scale. In step S1, the construction of the complex pore structure high-angle fracture type rock physical model comprises: s11. The equivalent modulus of the solid matrix is calculated The VRH average theory is used to calculate the equivalent bulk modulus of the solid matrix and the equivalent shear modulus of the solid matrix The expressions are as follows: ； ； wherein, is the volume fraction of the inclusion phase, is the number of inclusion phase classes, is the volume fraction of the inclusion phase, is the bulk modulus of the inclusion, is the shear modulus of the inclusion, is the equivalent bulk modulus of the solid matrix, is the equivalent shear modulus of the solid matrix; s12. The matrix dry skeleton containing multiple pores is constructed Using the K-T model to include the micro-pores into the solid matrix, for a solid matrix dry skeleton including randomly distributed micro-pores, the equivalent bulk modulus and shear modulus are expressed as follows: ； ； ； ； wherein, porosity of the micropores, and geometric factor of the micropores, equivalent bulk modulus of the solid matrix dry skeleton comprising randomly distributed micropores, equivalent shear modulus of the solid matrix dry skeleton comprising randomly distributed micropores; The differential effective medium theory is used to add different shapes of inter-particle pores to a solid matrix dry skeleton containing randomly distributed micro-pores to construct a matrix dry skeleton containing multiple pores, and the equivalent bulk modulus and shear modulus The expression is as follows: ； )； wherein, the porosity of the interparticulate porosity, and the geometric factor of the interparticulate porosity, the dry matrix skeleton equivalent bulk modulus of the matrix containing a plurality of pores, the dry matrix skeleton equivalent shear modulus of the matrix containing a plurality of pores; s13. The dry skeleton with micro-fracture development is constructed The micro-fractures in the directional arrangement are added to the matrix dry skeleton containing multiple pores by using the Hudson model to obtain a dry skeleton with micro-fracture development in the transverse isotropy, and the equivalent stiffness matrix is as follows: ； wherein, is the isotropic stiffness matrix of the matrix dry skeleton containing a plurality of pores, is a first order microcrack perturbation term, is a second order microcrack perturbation term; For convenience of model calculation, only the first-order micro-fracture perturbation term given by the Hudson model is calculated, and the specific expression is as follows: ； ； ； ； ； wherein, and is the Lame constant of the matrix dry skeleton containing a plurality of pores, is the fracture density of microfractures, and is a parameter depending on the fracture aspect ratio and the modulus of the background medium; s14. The dry rock skeleton containing fractures is constructed By using the Schoenberg linear sliding theory, two groups of macro high-angle vertical fractures are introduced on the basis of the dry skeleton with micro-fracture development, and the two groups of fractures are orthogonal in the azimuth, to construct a dry rock skeleton containing fractures; the equivalent stiffness matrix of the dry rock skeleton containing fractures is obtained through the stiffness matrix of the background medium and the fracture compliance matrix, the model assumes that the normal direction of the first group of fractures is along the x-axis, and the normal direction of the second group of fractures is along the y-axis, and the normal and tangential compliance matrix expressions are as follows: ； ； ； ； wherein, is the normal compliance of the first set of cracks, is the tangential compliance of the first set of cracks, is the normal compliance of the second set of cracks, is the tangential compliance of the second set of cracks, is the crack density of the first set of cracks, is the crack density of the second set of cracks, is the aspect ratio of the cracks of the first set of cracks, is the aspect ratio of the cracks of the second set of cracks, is a Poisson's ratio function; The equivalent stiffness matrix expression of the dry rock skeleton containing fractures is as follows: ； wherein is the crack compliance matrix, is the matrix composed of , , , ​ s15. The complex pore structure high-angle fracture type rock physical model is constructed The frequency dispersion and attenuation effects under the partial saturation condition are simulated by introducing the microscale and mesoscale double fluid flow mechanisms, the partially saturated fluid mixture is filled into the dry rock skeleton containing fractures to construct a fluid partially saturated tight sandstone rock, that is, a complex pore structure high-angle fracture type rock physical model, and the equivalent modulus expression of the model is as follows: ； ； ； wherein, is the frequency dependent partially saturated tight sandstone rock bulk modulus, is the dry rock matrix bulk modulus with fractures, is the frequency dependent partially saturated tight sandstone rock shear modulus, is the rock porosity, is the microfracture porosity, is the first set of fracture porosities, is the second set of fracture porosities, is the angular frequency, is the bulk modulus of the mixed fluid, is the imaginary unit, is the characteristic frequency determined by the fluid diffusion coefficient; The fluid partially saturated tight sandstone rock velocity under the theory is calculated, that is, the longitudinal wave velocity and the transverse wave velocity of the model, and the expressions are as follows: ； ； ； wherein, is the model compressional wave velocity, is the model shear wave velocity, is the solid matrix density, is the density of water, is the density of gas, is the water saturation, is the gas saturation.

3. The method of claim 2, wherein the method further comprises: In step s1, the inversion obtains two groups of fracture density parameters in the well , The step of the application specifically comprises: matching the actual shear wave velocity curve and the longitudinal wave velocity curve in the logging data with the theoretical shear wave velocity curve and the longitudinal wave velocity curve calculated by the model, constructing an objective function of the inversion, and defining the objective function as the error norm between the model prediction value and the measured value, and the expression is as follows: ； where, is the model parameter vector to be inverted, is the measured P-wave velocity in the well, is the measured S-wave velocity in the well, and is the weight coefficient, is the regularization term for introducing the prior information.

Images