Pre-stack inversion method for tight sand gas reservoirs with multiple sets of fracture density identification factors

Abstract

The application discloses a tight sandstone gas reservoir prestack inversion method of multiple sets of fracture density identification factors, and relates to the field of oil and gas geophysical exploration. The method firstly constructs a multi-scale pore-fracture structure partially saturated rock physical model to accurately describe the elastic response characteristics of tight sandstone containing high-angle multiple sets of fractures. Secondly, the method uses logging data to carry out multiple sets of fracture density rock physical inversion to realize quantitative evaluation of fracture density based on conventional logging. Thirdly, the method constructs two sets of fracture density identification factors based on Poisson impedance to effectively represent the development degree of specific fracture sets. Then, the method constructs AVO equation containing multiple sets of fracture density identification factors and solves the equation by using the improved grey wolf algorithm combined with Gaussian search threshold constraint to form a complete multiple sets of fracture density identification factor tight sandstone gas reservoir prestack inversion method. Practical data application shows that the method can accurately predict the development degree and spatial distribution of specific fracture sets of the reservoir, and provides technical support for tight sandstone gas exploration and development.

Description

Technical Field

[0001] This invention relates to the field of oil and gas geophysical exploration technology, and to a pre-stack inversion method for tight sandstone gas reservoirs using multiple sets of fracture density identification factors. Background Technology

[0002] Accurate description of fracture development characteristics is one of the core issues in the exploration and development of tight sandstone gas. High-angle fractures, a common type of fracture in tight sandstone gas reservoirs, are usually related to tectonic compression, regional faulting activity, or local stress fields. High-angle fractures not only provide effective reservoir space, but more importantly, they can significantly improve reservoir permeability, effectively connecting isolated pores and forming high-yield seepage channels, thus having a decisive impact on the development benefits of gas reservoirs. However, due to the strong heterogeneity of tight sandstone reservoirs and the multi-scale characteristics of fracture distribution, their geophysical response is weak, and their physical properties differ complexly from those of the surrounding rocks. Conventional seismic prediction methods are insufficient to quantitatively characterize their development density, spatial distribution, and fluid filling characteristics.

[0003] Currently, fracture reservoir prediction methods are mainly divided into two categories: post-stack attribute analysis and pre-stack anisotropic inversion, but both have obvious technical limitations. Specifically, these include: (1) Limitations of post-stack seismic attribute analysis methods. Traditional post-stack seismic attributes (such as curvature attributes, coherence attributes, variance volumes, etc.) mainly use the geometric discontinuity of seismic phase axes to indirectly infer fracture development zones. These methods have a certain ability to identify large-scale faults or fracture zones, but they have the following fundamental defects: they are not sensitive to small-scale fractures smaller than the seismic tuning thickness and cannot quantitatively evaluate fracture density; they cannot distinguish the directionality of fractures, especially the response to high-angle fractures is weak; they cannot provide information on the fluid properties of fractures and it is difficult to distinguish between gas-bearing fractures and water-bearing fractures; the attribute interpretation is ambiguous and lacks constraints from rock physics mechanisms. (2) Limitations of pre-stack anisotropic inversion methods. Fracture prediction methods based on amplitude variation with offset and azimuth (AVOA) utilize the anisotropic effect induced by fractures, and theoretically, fracture density and azimuth can be inverted simultaneously. However, existing methods suffer from the following technical bottlenecks: Conventional inversion methods are mostly based on the assumption of isotropic media or simple VTI (vertical-lateral isotropic) assumptions, making it difficult to accurately describe the complex anisotropic characteristics of reservoirs when multiple sets of fractures (especially high-angle fractures) coexist. Existing fracture petrophysics models (such as the Hudson model and the Schoenberg model) each have their own assumptions and applicable ranges, and a single model often cannot simultaneously consider the multi-scale characteristics of fractures, pore structure, and the influence of partially saturated fluids. Fracture density is a key parameter describing the degree of fracture development, but conventional logging data lacks a means to directly measure fracture density. How to establish a quantitative relationship between conventional logging data (P-wave velocity, S-wave velocity, density, porosity, clay content, etc.) and the density of multiple sets of fractures is a challenge in petrophysics research. Traditional elastic parameters (such as P-wave impedance, S-wave impedance, Lamé constant, etc.) have multiple solutions to the response of fracture fluids and are not sensitive enough to changes in fracture density, making it difficult to effectively distinguish between "gas-bearing high-angle fractures" and "non-fractured gas-bearing reservoirs" or "water-bearing fractures." Pre-stack seismic inversion is a complex optimization problem characterized by nonlinearity, high dimensionality, and multiple extrema. Traditional linearization methods (such as least squares) heavily rely on the initial model and are prone to getting trapped in local optima. While intelligent optimization algorithms (such as genetic algorithms, particle swarm optimization, and the standard Grey Wolf Optimization algorithm, GWO) have strong global search capabilities, they generally suffer from slow convergence, oscillations in later stages, and insufficient solution accuracy when handling high-dimensional, multi-parameter inversions. They lack effective local fine-grained search mechanisms and struggle to meet the efficiency and accuracy requirements of multi-channel seismic data inversion in practical production. Summary of the Invention

[0004] To address the technical challenges in existing technologies, such as the difficulty in quantitatively evaluating the development degree of multiple sets of high-angle fractures, low accuracy of rock physics modeling, and strong ambiguity in seismic inversion, this invention discloses a pre-stack inversion method for tight sandstone gas reservoirs using multiple sets of fracture density identification factors.

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

[0006] A pre-stack inversion method for tight sandstone gas reservoirs using multiple sets of fracture density identification factors, comprising:

[0007] Step 1: Data Acquisition and Preprocessing

[0008] We collected prestack seismic gather data, well logging data, mineral composition information, and geological stratification data for the study area; we performed amplitude preservation denoising, multiple removal, and dynamic correction preprocessing on the prestack seismic gather data to obtain P-wave amplitude data at different incident angles; and we performed environmental correction and depth matching on the well logging data.

[0009] Step 2: Partially saturated rock physics modeling of multi-scale porous and fractured structures.

[0010] A partially saturated rock physics model with high-angle multi-scale pore structure was constructed by comprehensively utilizing the Voigt-Reuss-Hill average formula, the improved Xu-White model, the differential equivalent medium theory (DEM), the Hudson model, the Schoenberg linear slip model, the improved White-Johnson partially saturated fluid model, and the Brown-Korringa anisotropic Gassmann equation.

[0011] Step 3: Multi-set fracture density rock physics inversion based on well logging data

[0012] Based on actual well logging data and a multi-scale pore structure partially saturated petrophysical model, multiple sets of nonlinear inversions of fracture density were constructed. The inversion process introduced a hybrid optimization algorithm that combines particle swarm optimization (PSO) with local gradient search.

[0013] Step 4: Constructing a crack density identification factor based on Poisson impedance

[0014] The concept of Poisson impedance is extended to fracture density identification. The variation characteristics of elastic parameters under different fracture density conditions of a specific fracture system are analyzed. Through parameter optimization driven by rock physics, an identification factor that is sensitive to the density of a specific fracture system but relatively insensitive to matrix lithology and pore fluid is constructed.

[0015] Step 5: Constructing the AVO approximation equation containing multiple sets of crack density identification factors.

[0016] The constructed crack density identification factor is directly embedded into the linear approximation of the Fatti high-precision AVO reflection coefficient, and the AVO approximation equation containing multiple sets of crack density identification factors is derived to establish a direct relationship between pre-stack seismic response and crack parameters.

[0017] Step 6: Pre-stack seismic inversion of tight sandstone gas reservoirs based on multiple sets of fracture density identification factors

[0018] Based on actual pre-stack seismic records and time-domain pre-stack synthetic seismic records, a least-norm 2 fitting is constructed, and a pre-stack seismic inversion objective function is defined. A Gaussian search threshold constraint-improved Grey Wolf Optimization Algorithm (GST-GWO) is proposed, which introduces adaptive Gaussian perturbation, an elite preservation strategy with threshold constraints, and a hybrid Cauchy-Gaussian mutation mechanism. The improved algorithm is embedded into the pre-stack inversion process to ensure that the inversion process can effectively escape local extrema caused by seismic noise and model errors, quickly and stably converge to the vicinity of the global optimum or a high-quality suboptimal solution, and achieve high-precision fitting.

[0019] Step 7: Comprehensive Interpretation of Crack Development Characteristics

[0020] A three-dimensional visualization analysis was performed on the multiple sets of fracture density identification factor data obtained by inversion. Combined with the results of geological structural interpretation, the development degree and spatial distribution pattern of multiple sets of fractures in tight sandstone gas reservoirs were quantitatively evaluated.

[0021] Furthermore, in step two, the process of partially saturated rock physics modeling of multi-scale porous structures is as follows:

[0022] 2.1 Construction of the mineral matrix of dense sandstone and calculation of its equivalent modulus

[0023] The equivalent bulk modulus of the mineral matrix was calculated using the Voigt-Reuss-Hill average. and equivalent shear modulus ,formula:

[0024] ;

[0025] ;

[0026] in, This is the equivalent bulk modulus of the mineral matrix; This is the equivalent shear modulus of the mineral matrix; This represents the volume fraction of clay minerals. This represents the bulk modulus of clay minerals. The shear modulus of clay minerals; This represents the volume fraction of quartz. The bulk modulus of quartz; This is the shear modulus of quartz;

[0027] 2.2 Construction of Porous Framework and Calculation of Equivalent Modulus in Tight Sandstone

[0028] The improved Xu-White model discretizes the pore space into flat, fracture-like pores and equiaxed pores to simulate intergranular pores and dissolution pores in the mineral matrix. The equivalent medium modulus of the porous framework in tight sandstone is calculated using the Kuster-Toksöz theory. The formula is:

[0029] ;

[0030] ;

[0031] ;

[0032] in, This is the equivalent bulk modulus of the mineral matrix; It is the equivalent bulk modulus of the porous skeleton; This is the equivalent shear modulus of the mineral matrix; Porosity refers to the porosity of crack-like pores; Geometric factors related to crack-like pores; Porosity refers to the porosity of equiaxed pores; Geometric factors related to equiaxed pores; This is the equivalent shear modulus of the mineral matrix; This is the equivalent shear modulus of the porous skeleton. The polarization factor associated with crack-like pores; The polarization factor associated with equiaxed pores; These are intermediate variables in the calculation process;

[0033] 2.3 Construction of the microfractured framework of dense sandstone and calculation of its equivalent modulus

[0034] The differential equivalent medium theory is used to incorporate randomly distributed microfractures into a porous framework to construct a dense sandstone framework containing microfractures. The formula for calculating the equivalent medium modulus is as follows:

[0035] ;

[0036] ;

[0037] in, The equivalent bulk modulus of the microcracked skeleton; This represents the volume fraction of microcracks. and A tensor characterizing the influence of microcrack shape; The equivalent bulk modulus of the microcrack; It is the equivalent bulk modulus of the porous skeleton; The equivalent shear modulus of the microcracked skeleton; The equivalent shear modulus of the microcrack; This is the equivalent shear modulus of the porous skeleton.

[0038] 2.4 Dry rock structure of dense sandstone and calculation of equivalent stiffness matrix

[0039] For multiple fracture systems, the influence of fractures on anisotropic elasticity is described by combining the Hudson model and the Schoenberg linear slip model. A set of high-angle fractures and a set of low-angle fractures are embedded in a dense sandstone skeleton with a microfracture framework as the background medium to construct a dry dense sandstone rock, in which the background medium is an isotropic medium.

[0040] For the Group of cracks, their crack density expression:

[0041] ;

[0042] in, Crack density; For the first Porosity of the crack group; For the first The average aspect ratio of the group of cracks;

[0043] The first group of cracks is defined as low-angle cracks, which are rotationally symmetric cracks arranged in a direction perpendicular to the X3 axis. The additional stiffness perturbation caused by the presence of low-angle cracks is given using the Hudson model. The expression for the first-order correction term:

[0044] ;

[0045] ;

[0046] in, The stiffness matrix of dry rock with low-angle cracks; The stiffness matrix of the background medium. The stiffness matrix is ​​the stiffness matrix for the first-order correction term of the added stiffness perturbation; The porosity of the first group of cracks; and Lamé constant of the background medium; , and This is the crack interaction term, used to characterize the crack state; , , , , and The Kronecker function; , , and The unit vector component perpendicular to the crack surface;

[0047] The second group of cracks is defined as high-angle cracks, which are rotationally symmetric cracks arranged in a direction perpendicular to the X1 axis. They are superimposed using the Schoenberg linear slip model. For the second group of cracks, its normal compliance... and tangential compliance The expression:

[0048] ;

[0049] ;

[0050] in, The crack density is for group 2. The equivalent bulk modulus of the microcracked skeleton; and The correlation coefficient for the crack infill material; The equivalent shear modulus of the microcracked skeleton;

[0051] The overall flexibility of the fractured dry rock is obtained by adding the flexibility tensor of the background medium to the additional flexibility tensors of multiple sets of fractures. Then, the inverse of this flexibility tensor is used to obtain the equivalent stiffness matrix of the dense sandstone dry rock. ,expression:

[0052] ;

[0053] in, Additional compliance for the second group of cracks;

[0054] 2.5 Characterization of Partially Saturated Fluid Effects

[0055] By incorporating fracture density parameters into the partially saturated model, an improved White-Johnson partially saturated fluid model is constructed. This model simplifies partially saturated rock into two types of periodic strata or spherical patches with different fluid saturations. The equivalent fluid bulk modulus considering the partially saturated fluid effect is calculated, expressed as:

[0056] ;

[0057] ;

[0058] ;

[0059] in, Gas saturation; The bulk modulus of the gas; The bulk modulus of water; The equivalent shear modulus of dry dense sandstone; It is a relaxation function; This is the Gassmann low-frequency limiting modulus; The relaxation strength; The imaginary unit; Angular frequency; The modified characteristic relaxation time is used to account for the effect of crack density. The characteristic relaxation time; The coupling coefficient is related to the crack direction; Crack density;

[0060] 2.6 Modeling and Calculation of Equivalent Stiffness Matrix for Partially Saturated Rock with Multi-Scale Pore and Fracture Structure in Tight Sandstone

[0061] Using the Brown-Korringa anisotropic Gassmann equations for fluid substitution, the anisotropic equivalent stiffness matrix of a partially saturated rock model with multi-scale porous structures is obtained, expressed as:

[0062] ;

[0063] in, The equivalent stiffness matrix of a partially saturated rock model with multi-scale porous and fractured structure in dense sandstone; The equivalent stiffness matrix of dense sandstone dry rock; and The equivalent stiffness matrix components of dense sandstone dry rock; This is the equivalent bulk modulus of the mineral matrix; and The Kronecker function; Porosity; The equivalent fluid bulk modulus that takes into account the effects of partially saturated fluid; This is the equivalent stiffness matrix of the mineral matrix.

[0064] Further, step three involves the process of multi-set fracture density rock physics inversion based on well logging data:

[0065] Using the model as the forward modeling operator and the P-wave and S-wave velocities from actual well logging data as constraints, the minimum fitting error between the measured well logging data and the model predictions is calculated. Then, multiple sets of fracture densities at each depth point are calculated, thereby inverting multiple sets of fracture density curves along the well depth. The inversion objective function is defined as the error between the measured well logging data and the model predictions. The expression is as follows:

[0066]

[0067] in, Forward modeling of the longitudinal wave velocity in rock physics; For logging actual P-wave velocity; The transverse wave velocity is simulated in the forward modeling of the rock physics model; For logging actual shear wave velocity; For regularization terms; The regularization coefficient is used. The crack density is for group 1. The crack density is for group 2.

[0068] A hybrid optimization algorithm combining particle swarm optimization and local gradient search is used for inversion solution. This algorithm combines the fast convergence characteristics of particle swarm optimization and the ability of local gradient search to escape local optima.

[0069] Furthermore, in step four, the process of constructing the crack density identification factor based on Poisson impedance is as follows:

[0070] Poisson impedance The general expression:

[0071] ;

[0072] ;

[0073] ;

[0074] in, For longitudinal wave impedance; Transverse wave impedance; The rotation factor; For longitudinal wave velocity; The transverse wave velocity; Density;

[0075] Based on a rock physics model, the elastic response characteristics of multiple sets of fractures under different fracture densities are analyzed, and the expression for the fracture density sensitivity function is defined:

[0076] ;

[0077] in, This refers to a certain elastic parameter under different crack densities in multiple sets of cracks; Crack density;

[0078] The optimal rotation coefficient is obtained by fitting the maximum correlation coefficient between the crack density sensitivity function and the Poisson impedance. ,expression:

[0079] ;

[0080] in, It is a Poisson impedance; It is a crack density sensitive function;

[0081] For two sets of cracks with different orientations, the concept of Poisson impedance is extended to crack density identification, defining two sets of crack density identification factors, expressed as:

[0082] ;

[0083] ;

[0084] in, The crack density identification factor for the first group of cracks; The optimal rotation coefficients for group 1; The crack density identification factor for the second group of cracks; The optimal rotation coefficients for group 2; For longitudinal wave impedance; Transverse wave impedance;

[0085] Step five involves the construction of the AVO approximation equation containing multiple sets of crack density identification factors:

[0086] The two sets of crack density identification factors based on Poisson impedance are extended to a differential expression, as shown in the following expression:

[0087] ;

[0088] ;

[0089] in, This represents the relative change in the crack density identification factor of the first group of cracks at the reflective interface. This represents the relative change in the crack density identification factor of the second group of cracks at the reflective interface. The relative rate of change of the crack density identification factor of the first group of cracks at the reflective interface; The relative rate of change of the crack density identification factor for the second group of cracks at the reflective interface; This represents the relative change in longitudinal wave impedance; This represents the relative change in transverse wave impedance; The optimal rotation coefficients for group 1; The optimal rotation coefficients for group 2; This represents the P-wave velocity; the current value is the same as the P-wave velocity measured in the well logging. same; This represents the shear wave velocity; the current value is the same as the shear wave velocity measured in the well logging. same;

[0090] Define intermediate variables , , , ,expression:

[0091] ;

[0092] ;

[0093] ;

[0094] ;

[0095] The simplified differential expression of the two sets of crack density identification factors based on Poisson impedance is as follows:

[0096] ;

[0097] ;

[0098] By solving the differential expressions for the two sets of crack density identification factors, the relative rate of change of the longitudinal wave impedance is obtained. The relative rate of change of transverse wave impedance ,expression:

[0099] ;

[0100] ;

[0101] Substituting the differential expressions of the two sets of crack density identification factors into the linear approximation of the Fatti high-precision AVO reflection coefficient, an AVO approximation equation containing multiple sets of crack density identification factors is derived. This yields the quantitative relationship between the PP wave reflection coefficient and the crack density identification factor, expressed as:

[0102] ;

[0103] in, The angle of incidence is denoted as .

[0104] Furthermore, in step six, during the direct inversion of pre-stack earthquakes, an objective function is constructed based on minimizing the fitting error between the actual pre-stack earthquake record and the synthetic earthquake record. The expression is:

[0105] ;

[0106] in, This is an actual earthquake record. For theoretically synthesized seismic records, N represents the number of incident angles.

[0107] The beneficial effects of this invention are that, compared with the prior art, its advantages are:

[0108] (1) The model is more consistent with geological reality

[0109] This invention constructs a multi-scale fracture partially saturated rock physics model, which overcomes the limitations of traditional single models. For the first time, it organically integrates background porosity (using the Xu-White and DEM models), multiple sets of high-angle fractures (using the Hudson and Schoenberg linear slip model), and partially saturated fluid effects (using an improved White-Johnson partially saturated fluid model). This model can more realistically reflect the microscopic pore structure and macroscopic fracture characteristics of tight sandstone reservoirs, laying a solid physical foundation for subsequent quantitative inversion.

[0110] (2) Quantification of parameter prediction

[0111] By using rock physics inversion, multiple sets of fracture densities can be quantitatively obtained using conventional well logging data, achieving a leap from "qualitative description" to "quantitative evaluation" and solving the problem of directly measuring fracture density in actual production. This provides crucial data support for extending fracture prediction from the well logging scale to three-dimensional space, and provides prior constraints and calibration basis for seismic inversion.

[0112] (3) The innovativeness and high sensitivity of the crack identification factor

[0113] The fracture density identification factor proposed in this invention, through parameter optimization driven by rock physics, achieves "purification" of information on specific fracture systems. It can distinguish fracture systems of different systems, making them sensitive to target fractures, while exhibiting good robustness to interference factors such as matrix lithology and pore fluids. This fundamentally improves the ability of seismic information to identify fractures and the reliability of fracture prediction.

[0114] (4) Solid theoretical foundation of inversion equation

[0115] The derived AVO approximation equation, containing multiple sets of crack density identification factors, establishes a direct link between pre-stack seismic response and crack parameters. Compared with traditional indirect prediction methods (first inverting elastic parameters and then calculating crack properties), this method avoids error propagation and accumulation, improving the accuracy and reliability of the inversion.

[0116] (5) The improved optimization algorithm has superior convergence performance.

[0117] The proposed hybrid optimization algorithm, combining Particle Swarm Optimization (PSO) with local gradient search, and the improved Grey Wolf Optimization (GST-GWO) algorithm with Gaussian search threshold constraints, offer two novel and improved intelligent solutions that balance global exploration and local exploitation capabilities. Compared to traditional algorithms, these algorithms offer faster convergence, higher solution accuracy, lower dependence on the initial model, and stronger resistance to seismic noise interference, making them highly suitable for handling complex seismic inversion problems.

[0118] (6) The technical process is systematic and highly practical.

[0119] This invention provides a complete technical system from well logging to seismic analysis, and from theory to application. Each step is interconnected, logically rigorous, highly operable, and easily integrated into existing seismic interpretation software platforms or self-developed systems, demonstrating promising prospects for industrial application.

[0120] (7) Excellent practical application results

[0121] Through practical data application verification, this technology can effectively identify the fracture properties of tight sandstone gas reservoirs. The planar distribution of fracture density has a good spatial matching relationship with the regional fracture system and structural curvature characteristics. High-density areas correspond to the vicinity of fracture zones and structural turning points, which is in line with the laws of geomechanics. The inversion results of fracture density identification factors have a good positive correlation with the actual gas production capacity. High fracture density areas correspond to high-yield well locations, providing a reliable basis for well location deployment and fracturing layer selection. Attached Figure Description

[0122] Figure 1 This invention provides a pre-stack inversion method for tight sandstone gas reservoirs based on multiple sets of fracture density identification factors.

[0123] Figure 2 This is a seismic two-way travel time (TWT) contour map of the target tight sandstone reservoir in the study area of ​​the application example of this invention. In the map, the pentagrams represent the locations of wells A, B, and C, the red pentagrams are gas wells, the black pentagrams are dry wells, and the black solid lines are seismic traces between wells.

[0124] Figure 3 This is an application example of the present invention, showing the results of rock physics inversion of multiple sets of fracture density based on well logging data in well A of the study area;

[0125] Figure 4 This is an application example of the present invention, showing the results of rock physics inversion of multiple sets of fracture density based on well logging data in well B of the study area.

[0126] Figure 5 This is the result of rock physics inversion of multiple sets of fracture density based on well logging data in well C of the study area, as an application example of this invention.

[0127] Figure 6The application example of this invention is the fracture density identification factor related parameter curve of well B in study area based on Poisson impedance;

[0128] Figure 7 This is an intersection diagram of the longitudinal wave impedance and the first group of fracture density identification factors in well B of the application example of this invention. The color scale represents the first group of fracture density.

[0129] Figure 8 This is an intersection diagram of the longitudinal wave impedance and the first group of fracture density identification factors in well B of the application example of this invention. The color scale represents the second group of fracture density.

[0130] Figure 9 This is a cross-plot of the longitudinal wave impedance and the second group of fracture density identification factors in well B of the application example of this invention. The color scale represents the second group of fracture density.

[0131] Figure 10 This is an intersection diagram of the longitudinal wave impedance and the second group of fracture density identification factors in well B of the application example of this invention. The color scale represents the first group of fracture density.

[0132] Figure 11 This is a well profile of the first group of fracture density identification factors in the application example study area of ​​this invention;

[0133] Figure 12 This is a well profile of the second group of fracture density identification factors in the application example study area of ​​this invention;

[0134] Figure 13 This is a planar diagram of the crack density identification factor of the first group in the target layer of the study area in the application example of this invention;

[0135] Figure 14 This is a plan view of the crack density identification factor of the second group in the target layer of the application example study area of ​​the present invention. Detailed Implementation

[0136] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of them. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0137] This invention discloses a pre-stack inversion method for tight sandstone gas reservoirs using multiple sets of fracture density identification factors. It establishes a multi-scale petrophysical model capable of describing the coupling effect between high-angle multiple sets of fractures and partially saturated fluids; develops a quantitative method to directly extract multiple sets of fracture density parameters from conventional logging data; constructs fracture identification factors sensitive to fracture density and capable of suppressing lithological interference; and develops a nonlinear optimization algorithm that enables direct inversion of reservoir parameters from pre-stack seismic data. Furthermore, it constructs an integrated, high-precision, and practical inversion technology and process based on petrophysical mechanisms, effectively integrating logging and seismic data to accurately quantify the development characteristics of high-angle multiple sets of fractures in tight sandstone gas reservoirs.

[0138] This invention uses a partially saturated petrophysical model of multi-scale pore and fracture structure in tight sandstone as its foundation, multiple sets of fracture density identification factors based on Poisson impedance as its bridge, the AVO approximation equation containing multiple sets of fracture density identification factors as its engine, and an improved intelligent optimization inversion algorithm as its driving force to construct an integrated inversion framework of "rock physics-logging-seismic" to achieve direct and quantitative prediction of high-angle multiple sets of fracture density parameters in tight sandstone gas reservoirs. This provides key data support for the prediction of "sweet spot" areas, optimization of horizontal well trajectories, and reserve calculation in tight sandstone gas reservoirs.

[0139] Figure 1 The flowchart below shows the pre-stack inversion method for tight sandstone gas reservoirs based on multiple sets of fracture density identification factors, which includes the following steps:

[0140] (I) Multi-scale porous structure partially saturated rock physical modeling

[0141] Based on the theory of multi-scale porosity and fracture coexistence, the pore space of tight sandstone reservoirs is divided into three scales: matrix porosity (nano-micron scale), mainly intergranular pores and dissolution pores; microfractures (micron-millimeter scale), usually randomly distributed, with variable occurrence, approximately isotropic, but significantly reducing the overall stiffness of the rock; and macrofractures (millimeters and above), with a clear dominant orientation (such as high-angle tectonic fractures). Among them, matrix porosity and microfractures constitute the background isotropic elastic framework of the rock and are the main reservoir space of the matrix, while macrofractures constitute the main seepage channels. The rock physics model of this invention adopts a step-by-step construction algorithm, treating the tight sandstone reservoir as a complex system composed of rock matrix, pore space, multiple sets of high-angle fractures, and partially saturated fluids.

[0142] 1.1 Construction of the mineral matrix and calculation of equivalent modulus of dense sandstone

[0143] Based on well logging data, the main mineral components of tight sandstone reservoirs are quartz and clay. The equivalent bulk modulus of the mineral matrix can be calculated using the volume fractions of these mineral components obtained from well logging interpretation, and the Voigt-Reuss-Hill average can be applied. and equivalent shear modulus The formula is as follows:

[0144] (1);

[0145] (2);

[0146] in, This is the equivalent bulk modulus of the mineral matrix; This is the equivalent shear modulus of the mineral matrix; This represents the volume fraction of clay minerals. This represents the bulk modulus of clay minerals. The shear modulus of clay minerals; This represents the volume fraction of quartz. The bulk modulus of quartz; This is the shear modulus of quartz.

[0147] 1.2 Construction of Porous Framework and Calculation of Equivalent Modulus in Tight Sandstone

[0148] The elastic properties of the porous framework of tight sandstone are controlled by porosity and pore aspect ratio. Tight sandstone has a complex pore structure. To address the influence of pore aspect ratio distribution on the elastic modulus, an improved Xu-White model is used. The improved Xu-White model discretizes the pore space into flat, fracture-like pores (high aspect ratio) and equiaxed pores (low aspect ratio) to simulate intergranular pores and dissolution pores in the mineral matrix. The equivalent medium modulus of the porous framework of tight sandstone is calculated using the Kuster-Toksöz theory, as shown in the following formula:

[0149] (3);

[0150] (4);

[0151] (5);

[0152] in, This is the equivalent bulk modulus of the mineral matrix; It is the equivalent bulk modulus of the porous skeleton; This is the equivalent shear modulus of the mineral matrix; Porosity refers to the porosity of crack-like pores; Geometric factors related to crack-like pores; Porosity refers to the porosity of equiaxed pores; Geometric factors related to equiaxed pores; This is the equivalent shear modulus of the mineral matrix; This is the equivalent shear modulus of the porous skeleton. The polarization factor associated with crack-like pores; The polarization factor associated with equiaxed pores; These are intermediate variables in the calculation process.

[0153] 1.3 Construction of the microfractured framework of dense sandstone and calculation of its equivalent modulus

[0154] Tight sandstone often contains isolated microfractures that are not connected to the main pore space and are usually randomly distributed in the background medium. Using the differential equivalent medium (DEM) theory, these randomly distributed microfractures are added to the porous framework, forming a dry "matrix-fractured pores-equiaxed pores-microfractures" equivalent medium, thus constructing a tight sandstone framework containing microfractures. The formula for calculating the equivalent medium modulus is as follows:

[0155] (6);

[0156] (7);

[0157] in, The equivalent bulk modulus of the microcracked skeleton; This represents the volume fraction of microcracks. and A tensor characterizing the influence of microcrack shape; The equivalent bulk modulus of the microcrack; It is the equivalent bulk modulus of the porous skeleton; The equivalent shear modulus of the microcracked skeleton; The equivalent shear modulus of the microcrack; It is the equivalent shear modulus of the porous skeleton.

[0158] 1.4 Dry rock structure of dense sandstone and calculation of equivalent stiffness matrix

[0159] In actual tight sandstone reservoirs, multiple sets of fracture systems with different orientations and high angles are often developed. High-angle fractures are not only channels for oil and gas migration but also the main reservoir spaces and seepage channels, and their development level directly determines the distribution of "sweet spots." To address the fracture development characteristics of tight sandstone gas reservoirs, this invention constructs a rock physical model containing two sets of fracture systems. To comprehensively analyze the fracture characteristics of actual reservoirs, a set of high-angle fracture systems and a set of low-angle fracture systems are introduced into the model. For multiple fracture systems, the influence of fractures on anisotropic elasticity is described by combining the Hudson model and the Schoenberg linear slip model. Using a tight sandstone microfractured framework as the background medium, a set of high-angle fractures and a set of low-angle fractures are embedded to construct a dry tight sandstone rock, where the background medium is isotropic.

[0160] For the Group of cracks, their crack density The expression is as follows:

[0161] (8);

[0162] in, Crack density; For the first Porosity of the crack group; For the first The average aspect ratio of the group of cracks;

[0163] The first group of cracks is defined as low-angle cracks, which are rotationally symmetric cracks arranged in a direction perpendicular to the X3 axis. The additional stiffness perturbation caused by the presence of low-angle cracks is given using the Hudson model. Its first-order correction term is expressed as follows:

[0164] (9);

[0165] (10);

[0166] in, The stiffness matrix of dry rock with low-angle cracks; The stiffness matrix of the background medium. The stiffness matrix is ​​the stiffness matrix for the first-order correction term of the added stiffness perturbation; The porosity of the first group of cracks; and Lamé constant of the background medium; , and This is the crack interaction term, used to characterize the crack state; , , , , and The Kronecker function; , , and The unit vector component perpendicular to the crack surface;

[0167] When the crack density is low, the Hudson first-order approximation can be used. When the crack density is high, in order to avoid the high-density instability of the Hudson model, the self-consistent approximation (SCA) model is introduced for correction, and the cracks and the background medium are regarded as equivalent media for iterative solution.

[0168] The second group of cracks is defined as high-angle cracks, which are rotationally symmetric cracks arranged in a direction perpendicular to the X1 axis, and are superimposed using the Schoenberg linear slip model. This model treats the cracks as interfaces of displacement discontinuity and uses the crack compliance tensor to describe their mechanical effects. For the second group of cracks, its normal compliance... and tangential compliance The expression is as follows:

[0169] (11);

[0170] (12);

[0171] in, The crack density is for group 2. The equivalent bulk modulus of the microcracked skeleton; and The correlation coefficient for the crack infill material; It represents the equivalent shear modulus of the microcracked skeleton.

[0172] The overall flexibility of the fractured dry rock is obtained by adding the flexibility tensor of the background medium to the additional flexibility tensors of multiple sets of fractures. The inverse of this flexibility tensor yields the equivalent stiffness matrix of the fractured dry rock, which is the equivalent stiffness matrix of the dense sandstone dry rock. The expression is as follows:

[0173] (13);

[0174] in, Additional flexibility for the second group of cracks.

[0175] 1.5 Characterization of Partially Saturated Fluid Effects

[0176] In tight sandstone gas reservoirs, both natural gas and formation water coexist within the pore space. Due to the heterogeneity of the pore structure, these two fluids are distributed in patchy patterns. This non-uniform saturation characteristic of gas and water is known as a non-uniform partial saturation state (patch saturation) at the mesoscale. When seismic waves propagate, they create a pressure gradient between the saturated patches and the surrounding partially saturated areas, leading to fluid flow (wave-induced flow), which in turn produces significant seismic wave velocity dispersion and energy attenuation.

[0177] For the non-uniform gas-water saturation characteristics of tight sandstone gas reservoirs, the dispersion and attenuation caused by wave-induced fluid flow can be described by the partial saturation theory at the mesoscale. In cases with multiple fractures, the presence of fractures significantly affects local fluid flow. This invention introduces the fracture density parameter into the partial saturation model, constructing an improved White-Johnson partially saturated fluid model. This model simplifies partially saturated rocks into two types of periodic strata or spherical patches with different fluid saturations, and calculates the equivalent fluid bulk modulus considering the partial saturation fluid effect, as expressed below:

[0178] (14);

[0179] (15);

[0180] (16);

[0181] in, Gas saturation; The bulk modulus of the gas; The bulk modulus of water; The equivalent shear modulus of dry dense sandstone; It is a relaxation function; This is the Gassmann low-frequency limiting modulus; The relaxation strength; The imaginary unit; Angular frequency; The modified characteristic relaxation time is used to account for the effect of crack density. The characteristic relaxation time; The coupling coefficient is related to the crack direction; The crack density is given.

[0182] 1.6 Modeling and Calculation of Equivalent Stiffness Matrix for Partially Saturated Rock with Multi-Scale Pore and Fracture Structure in Tight Sandstone

[0183] By combining dry, compact sandstone with partially saturated fluid, a multi-scale porous, partially saturated rock model of compact sandstone was constructed. This invention utilizes the Brown-Korringa anisotropic Gassmann equation for fluid substitution to obtain the anisotropic equivalent stiffness matrix of the multi-scale porous, partially saturated rock model, expressed as:

[0184] (17);

[0185] in, The equivalent stiffness matrix of a partially saturated rock model with multi-scale porous and fractured structure in dense sandstone; The equivalent stiffness matrix of dense sandstone dry rock; and The equivalent stiffness matrix components of dense sandstone dry rock; This is the equivalent bulk modulus of the mineral matrix; and The Kronecker function; Porosity; The equivalent fluid bulk modulus that takes into account the effects of partially saturated fluid; This is the equivalent stiffness matrix of the mineral matrix.

[0186] This model establishes a complete mapping relationship from microscopic pore structure (pore type, porosity, pore aspect ratio), crack characteristics (crack system, density, attitude) to macroscopic elastic response (equivalent anisotropic stiffness matrix), laying the physical foundation for subsequent crack density inversion.

[0187] (II) Rock Physical Inversion of Multiple Fracture Density Based on Well Logging Data

[0188] The idea of ​​multi-set fracture density rock physics inversion based on well logging data is to construct a mapping relationship from fracture density to P-wave and S-wave velocities, with a rock physics model as the core and conventional well logging data as the main observation data.

[0189] This invention is based on a partially saturated petrophysical model of multi-scale pore and fracture structure in tight sandstone. Using the model as a forward modeling operator and the P-wave and S-wave velocities from actual well logging data as constraints, it calculates the minimum fitting error between the measured well logging data and the model's predicted values. This allows for the calculation of multiple sets of fracture densities at each depth point, thereby inverting multiple sets of fracture density curves along the well depth. The objective function for inversion is defined as the error between the measured well logging data and the model's predicted values. The expression is as follows: ;

[0190] in, Forward modeling of the longitudinal wave velocity in rock physics; For logging actual P-wave velocity; The transverse wave velocity is simulated in the forward modeling of the rock physics model; For logging actual shear wave velocity; For regularization terms; The regularization coefficient is used. The crack density is for group 1. The crack density is for group 2.

[0191] Considering the highly nonlinear nature of the rock physics model and the fact that solving the objective function is a nonlinear optimization problem, this invention employs a hybrid optimization algorithm combining Particle Swarm Optimization (PSO) and Local Gradient Search (LTS) for inversion. This algorithm combines the fast convergence of PSO with the ability of LGS to escape local optima. During iteration, it accepts a degraded solution with a certain probability to avoid getting trapped in local minima. The main algorithm flow is as follows: initialize the particle swarm (initialize model parameters); calculate the forward response and objective function value of each particle; update the individual optimum and global optimum; iterate until the convergence condition is met; use the global optimum as the initial value and perform LGS to improve accuracy.

[0192] Through the above inversion, multiple sets of fracture densities were successfully quantitatively inverted using conventional well logging data, obtaining fracture density values ​​at each depth point of the target well location, and thus obtaining multiple sets of continuous fracture density curves, providing prior constraints and calibration data for subsequent seismic inversion.

[0193] (III) Construction of Crack Density Identification Factor Based on Poisson Impedance

[0194] Poisson impedance, a fluid identification sensitive parameter developed in recent years, is defined as a linear combination of P-wave impedance and S-wave impedance, achieving coordinate rotation transformation to enhance sensitivity to changes in specific rock physical properties. This invention extends the concept of Poisson impedance to the field of fracture identification, aiming to construct a fracture density identification factor that is sensitive to fracture density and relatively stable to changes in lithology and fluids, serving as a "characteristic parameter" connecting rock physical parameters and seismic reflection.

[0195] Poisson impedance The general expression is as follows:

[0196] (19);

[0197] (20);

[0198] (twenty one);

[0199] in, For longitudinal wave impedance; Transverse wave impedance; The rotation factor; For longitudinal wave velocity; The transverse wave velocity; Density;

[0200] Based on a rock physics model, the elastic response characteristics of multiple sets of fractures under different fracture densities are analyzed. A fracture density sensitivity function is defined, expressed as follows:

[0201] (twenty two);

[0202] in, For multiple sets of cracks with different crack densities, a certain elastic parameter is used, such as longitudinal wave impedance, transverse wave impedance, longitudinal wave to transverse wave velocity ratio, Poisson's ratio, anisotropy parameter, etc. The value represents the fracture density. Analysis shows that Poisson's ratio and the P-wave / S-wave velocity ratio are relatively sensitive to fracture density, while the P-wave / S-wave impedance ratio is significantly affected by lithology.

[0203] The optimal rotation coefficient is obtained by fitting the maximum correlation coefficient between the crack density sensitivity function and the Poisson impedance. The expression is as follows:

[0204] (twenty three);

[0205] in, It is a Poisson impedance; This is a crack density sensitivity function.

[0206] For two sets of cracks with different orientations, the concept of Poisson impedance is extended to crack density identification, and two sets of crack density identification factors are defined, as follows:

[0207] (twenty four);

[0208] (25);

[0209] in, The crack density identification factor for the first group of cracks; The optimal rotation coefficients for group 1; The crack density identification factor for the second group of cracks; The optimal rotation coefficients for group 2; For longitudinal wave impedance; This is the transverse wave impedance.

[0210] (iv) Construction of the AVO approximation equation containing multiple sets of crack density identification factors

[0211] To achieve direct inversion of fracture parameters based on pre-stack seismic data, it is necessary to establish an AVO approximation equation containing a fracture density identification factor. Traditional AVO equations (such as the Zoeppritz equation, the Fatti reflection coefficient approximation, and the Aki-Richards reflection coefficient approximation) only include background elastic parameters and cannot directly establish a relationship with fracture density. This invention introduces a fracture density identification factor based on Poisson impedance into the AVO equation, constructing an AVO approximation equation containing multiple sets of fracture density identification factors. This achieves a direct correlation between reservoir parameters and seismic response. Based on anisotropic elasticity theory, an approximate formula for the PP wave reflection coefficient containing the fracture density identification factor is then derived.

[0212] The two sets of crack density identification factors based on Poisson impedance are extended to a differential expression, as shown below:

[0213] (26);

[0214] (27);

[0215] in, This represents the relative change in the crack density identification factor of the first group of cracks at the reflective interface. This represents the relative change in the crack density identification factor of the second group of cracks at the reflective interface. The relative rate of change of the crack density identification factor of the first group of cracks at the reflective interface; The relative rate of change of the crack density identification factor for the second group of cracks at the reflective interface; This represents the relative change in longitudinal wave impedance; This represents the relative change in transverse wave impedance; The optimal rotation coefficients for group 1; The optimal rotation coefficients for group 2; This represents the P-wave velocity; the current value is the same as the P-wave velocity measured in the well logging. same; This represents the shear wave velocity; the current value is the same as the shear wave velocity measured in the well logging. same.

[0216] Because the derivation involves many parameters, intermediate variables are defined to simplify the formula derivation process. , , , The expression is as follows:

[0217] (28);

[0218] (29);

[0219] (30);

[0220] (31);

[0221] The simplified differential expressions for the two sets of crack density identification factors based on Poisson impedance are as follows:

[0222] (32);

[0223] (33);

[0224] By solving the differential expressions for the two sets of crack density identification factors, the relative rate of change of the longitudinal wave impedance is obtained. The relative rate of change of transverse wave impedance The expression is as follows:

[0225] (34);

[0226] (35);

[0227] Substituting the differential expressions of the two sets of crack density identification factors into the linear approximation of the Fatti high-precision AVO reflection coefficient, an AVO approximation equation containing multiple sets of crack density identification factors is derived. The quantitative relationship between the PP wave reflection coefficient and the crack density identification factor can then be obtained, as shown in the following expression:

[0228] (36);

[0229] in, The angle of incidence is denoted as .

[0230] This equation establishes an explicit relationship between pre-stack seismic gathers at different incident angles and two sets of fracture density identification factors, clarifying the relationship between the PP wave reflection coefficient and the two sets of fracture density identification factors. , and density changes The quantitative relationship between them provides a theoretical basis for the subsequent direct inversion of crack parameters.

[0231] (V) Pre-stack inversion method for tight sandstone gas reservoirs based on multiple sets of fracture density identification factors

[0232] Based on the AVO approximation equation containing multiple sets of crack density identification factors, the reflection coefficient of PP wave under different incident angles is calculated. By convolving it with the seismic wavelet spectrum, the pre-stack synthetic seismic record in the frequency domain can be obtained. Then, the pre-stack synthetic seismic record in the time domain can be obtained by using the inverse Fourier transform.

[0233] In pre-stack earthquake direct inversion, the objective function is constructed based on minimizing the error between the actual pre-stack earthquake record and the synthetic earthquake record, as shown in the following expression:

[0234] (37);

[0235] in, This is an actual earthquake record. For theoretically synthesized seismic records, N represents the number of incident angles.

[0236] Pre-stack seismic inversion of tight sandstone gas reservoirs based on multiple sets of fracture density identification factors is a highly nonlinear, multi-parameter, and multi-extremum optimization problem. Optimization algorithms can be introduced to improve the efficiency of optimizing complex objective functions. This invention proposes a Gaussian search threshold constraint-improved Grey Wolf Optimization Algorithm (GST-GWO) to solve the problem of direct reservoir parameter inversion. The traditional Grey Wolf Optimization Algorithm (GWO) simulates the social hierarchy and hunting behavior of grey wolf packs, possessing advantages such as simple structure, few parameters, and fast convergence speed. However, it suffers from drawbacks such as being prone to getting trapped in local optima and insufficient convergence accuracy in later stages. To improve the global search capability and convergence accuracy of the GWO algorithm, this invention introduces a Gaussian search threshold constraint mechanism. This mechanism includes three core improvements:

[0237] (1) Adaptive Gaussian perturbation: In the process of position update, an adaptive Gaussian perturbation term is introduced to enhance population diversity and avoid premature convergence. This design enables the algorithm to have strong exploration capabilities in the early stage and gradually focus on local fine search in the later stage.

[0238] (2) Threshold-constrained elite retention algorithm: Define convergence state threshold and diversity threshold. When the population convergence (measured by the average distance between individuals in the population and α wolves) is lower than 0 and the population diversity is lower than 0, the threshold constraint mechanism is triggered. This algorithm reactivates the search capability when the population is trapped in a local extreme value.

[0239] (3) Hybrid Cauchy-Gaussian Mutation: In the later stages of iteration, a hybrid mutation operation is applied to the current best individual to increase the probability of escaping the local optimum. The thick-tailed characteristic of Cauchy mutation helps to generate large-step perturbations, while Gaussian mutation provides local fine-grained search capabilities. The process of solving the pre-stack seismic inversion problem using the improved Gray Wolf Optimization Algorithm (GST-GWO) with Gaussian search threshold constraint is as follows: Initialize the gray wolf population (randomly generate candidate model parameters); calculate the fitness (objective function value) of each individual; determine α, β, and δ wolves (the three individuals with the best fitness); update the population according to the improved position update formula; apply the Gaussian search threshold constraint mechanism to determine whether a perturbation is triggered; if the termination condition is met (reaching the maximum number of iterations or the fitness change is less than the threshold), output α wolf as the optimal solution, i.e., the optimal individual position; post-process the inversion results and output two sets of fracture density identification factors.

[0240] This method has the following characteristics:

[0241] (1) Provides a multi-scale porous structure partially saturated rock physics model: integrates the Voigt-Reuss-Hill average formula, the improved Xu-White model, the differential equivalent medium theory (DEM), the Hudson model, the Schoenberg linear slip model, the improved White-Johnson partially saturated fluid model and the Brown-Korringa anisotropic Gassmann equation to construct a multi-scale porous structure partially saturated rock physics model with high angles. It builds a unified rock physics theoretical framework that can simultaneously describe multiple sets of high-angle high-fractures, matrix porosity (intergranular pores and dissolution pores), microfractures and gas-water non-uniform saturation state, and solves the problem that existing models are insufficient in describing the multi-scale porous structure-fluid coupling effect.

[0242] (2) A multi-set fracture density rock physics inversion based on well logging data is provided: a multi-set fracture density nonlinear inversion is constructed based on actual well logging data and a multi-scale pore structure partially saturated rock physics model. The inversion process introduces a hybrid optimization algorithm that combines particle swarm optimization (PSO) and local gradient search to achieve quantitative prediction of multi-set fracture density, providing prior constraints and calibration basis for seismic inversion.

[0243] (3) A fracture density identification factor based on Poisson impedance is proposed: For two sets of fracture systems (one of which is high-angle fractures), the concept of Poisson impedance is extended to fracture density identification. The elastic parameter variation characteristics of a specific fracture system under different fracture density conditions are analyzed. Through parameter optimization driven by rock physics, an identification factor that is sensitive to the density of a specific fracture system but relatively insensitive to matrix lithology and pore fluids is constructed. This achieves the "purification" of information on a specific fracture system, enabling the differentiation of fracture systems of different systems, making it sensitive to the target fracture, while exhibiting good robustness to interference factors such as matrix lithology and pore fluids. This fundamentally improves the ability of seismic information to identify fractures and the reliability of fracture prediction.

[0244] (4) A method is proposed to construct an AVO approximation equation containing multiple sets of fracture density identification factors: the constructed fracture density identification factors are directly embedded into the linear approximation of the Fatti high-precision AVO reflection coefficient, and the derived AVO approximation equation containing multiple sets of fracture density identification factors establishes a direct relationship between pre-stack seismic response and fracture parameters. Compared with traditional indirect prediction methods, this method avoids error propagation and accumulation, improves the accuracy and reliability of inversion, and provides a forward modeling basis for the direct inversion of pre-stack seismic data.

[0245] (5) Providing an improvement to traditional optimization algorithms: While traditional particle swarm optimization (PSO) possesses global search capabilities, it suffers from slow convergence, late-stage oscillations, and premature convergence. This is addressed by combining PSO with local gradient search to construct a hybrid optimization algorithm that integrates PSO and local gradient search for inversion solutions. This algorithm combines the fast convergence of PSO with the ability of local gradient search to escape local optima. During iteration, it accepts degraded solutions with a certain probability to avoid getting trapped in local minima. To address the premature convergence and insufficient global search capabilities of traditional gray wolf optimization (GGSO), a Gaussian search threshold constraint mechanism is introduced to improve the convergence speed and solution accuracy of nonlinear inversion, achieving stable, efficient, and direct inversion of reservoir parameters.

[0246] (6) Provide a complete technical process for pre-stack inversion of tight sandstone gas reservoirs with multiple sets of fracture density identification factors: verify its effectiveness and superiority through actual data application, perform three-dimensional visualization analysis on the data volume of multiple sets of fracture density identification factors obtained by inversion, and quantitatively evaluate the development degree and spatial distribution pattern of multiple sets of fractures in tight sandstone gas reservoirs in combination with the results of geological structural interpretation.

[0247] Application examples

[0248] The study area is located on the eastern edge of the Ordos Basin and in the northwestern part of the Shanxi-Western Shanxi flexural belt. The overall thickness of the sedimentary strata within the basin is 4000–6000 m, and hydrocarbon resources are mainly hosted in the Upper Paleozoic strata. The stratigraphy in this area generally trends northeast to southwest, and the tectonic evolution process is consistent with the overall trend of the basin. As a typical tectonic-sedimentary transition zone, this area is of great research significance in terms of fracture development, sedimentary response, and hydrocarbon accumulation.

[0249] Figure 2 This is a seismic two-way travel time (TWT) contour map of the target tight sandstone reservoir in the study area, covering a time range of 920 ms to 980 ms. The study area comprises 200 inline lines and 700 xline lines, with an inline spacing of 40 m and an xline spacing of 20 m, covering a total area of ​​approximately 120 km². 2 The seismic two-way travel time (TWT) contour map of the target tight sandstone reservoir in the study area shows that the overall geological structure in the study area is relatively gentle. The location of three wells (Well A, Well B, and Well C) is marked with a pentagram in the map. Wells A and C are dry wells, and Well B is a gas producing well. The black solid line represents the seismic tunnel between the wells.

[0250] A multi-set fracture density petrophysical inversion method based on well logging data was used to quantitatively identify two sets of fractures at well locations in the study area using well logging data from three wells (A, B, and C). For each sampling point in the well, well logging data and mineral composition information were used as input data. Based on a multi-scale partially saturated petrophysical model of pore and fracture structure, the preset fracture density values ​​at each depth point and the corresponding P-wave and S-wave velocities simulated by the model forward model were calculated. The simulated velocities were fitted with the measured velocities. When the fitting difference between the two reached its maximum, the two sets of optimal fracture density parameters for that depth point were obtained. This allowed for the inversion of the two sets of fracture densities and their corresponding P-wave and S-wave velocity curves at the well logging scale. During the inversion process, a hybrid optimization algorithm combining particle swarm optimization (PSO) and local gradient search was introduced for the inversion solution.

[0251] The results of multiple sets of fracture density petrophysical inversion based on well logging data in well AC in the study area are as follows: Figure 3-5 As shown in the figure, the logging data of well AC includes gamma rays, P-wave velocity, S-wave velocity, density, porosity, water saturation, and mineral composition information. It also presents two sets of fracture densities and their corresponding P-wave and S-wave velocities obtained based on the method of this invention. Well B is a gas well, and its logging data exhibits a "three lows and one high" characteristic: low gamma rays, low density, low water saturation, and high velocity, indicating that the reservoir in the study area has good geological response characteristics of tight sandstone gas reservoirs. The results show that the P-wave velocities obtained from the inversion in the three wells have a high degree of fit with the P-wave velocities measured in the wells, and the simulated S-wave velocities also show a high degree of consistency with the S-wave velocities measured in the wells. This verifies the excellent applicability of the multi-scale pore and fracture structure partially saturated petrophysical model proposed in this invention to the target tight sandstone reservoir in the study area, and corroborates the accuracy of the two sets of fracture densities obtained from the inversion. The two sets of fracture density curves exhibit distinct zonation characteristics at different depths, which are consistent with the changes in lithology, porosity, and fluid-bearing properties. This not only verifies the accuracy and applicability of the method of this invention in the quantitative prediction of fracture density in tight sandstone reservoirs, but also provides reliable technical support for the characterization and development of fractured reservoir parameters.

[0252] Based on the measured data from Well B in the study area and the above two sets of fracture density inversion results, the elastic parameter variation characteristics of a specific fracture system under different fracture density conditions were analyzed through parameter optimization driven by rock physics. Based on the Poisson impedance equation, the optimal rotation coefficient was calculated, and an identification factor that is sensitive to the density of a specific fracture system but relatively insensitive to matrix lithology and pore fluid was constructed. Figure 6The curves of the fracture density identification factor based on Poisson impedance for well B in the study area are shown. The results indicate that, for the target well location, the fracture density identification factor of a specific fracture group shows a significant positive correlation with the corresponding fracture density; that is, high fracture density areas correspond to high fracture density identification factor values. The correlation between the fracture density identification factor and fracture density is poor for different fracture groups, enabling the identification of fracture systems from different groups. This verifies the effectiveness of the fracture density identification factor in characterizing fracture development in tight reservoirs.

[0253] Figure 7-10 A cross-plot of longitudinal wave impedance and two sets of crack density identification factors is presented, with color scales indicating the two sets of crack densities. , The results indicate that the crack density identification factor of a specific crack group is strongly correlated with the longitudinal wave impedance, regardless of which crack density is represented by the color code. Figure 7 and Figure 9 The results show that the crack density identification factor of a specific crack system has a strong indicative effect on the corresponding crack density. Different crack density regions show obvious separation, and high crack density (orange) corresponds to the region with high crack density identification factor value, which further verifies the superiority of the crack density identification factor parameter in identifying the degree of crack development. Figure 8 and Figure 10 The results show that the fracture density identification factor for a specific fracture system has a low correlation with the density of non-corresponding fractures and weak distinguishability between different fracture densities. This verifies that the fracture density identification factor based on Poisson impedance constructed in this invention is only sensitive to the density of the corresponding fracture system and is insensitive to the density of non-corresponding fractures. In summary, the fracture density identification factor based on Poisson impedance constructed in this invention achieves "purification" of information on specific fracture systems, can distinguish fracture systems of different systems, making it sensitive to the target fracture, and has good robustness to interference factors such as matrix lithology and pore fluids. This fundamentally improves the ability of seismic information to identify fractures and can effectively characterize the degree of fracture development. It can be used as a parameter to characterize the fracture development characteristics of tight sandstone gas reservoirs.

[0254] The constructed Poisson impedance-based fracture density identification factor is embedded into the linear approximation of the Fatti high-precision AVO reflection coefficient, deriving an AVO approximation equation containing multiple sets of fracture density identification factors. This leads to the PP wave reflection coefficient formula, which is then compared with pre-stack synthetic seismic records generated in the time domain of a 25Hz Ricker wavelet. Input data consists of pre-stack seismic records with incident angles ranging from 5° to 30° in the study area. The upper and lower limits of the fracture parameter search space are set to the left and right 50% of the prior information values. A pre-stack inversion method for tight sandstone gas reservoirs based on multiple sets of fracture density identification factors is used to invert two sets of fracture density identification factor data volumes. During optimization, a dynamic weighting and local search algorithm is employed, and an adaptive convergence factor is introduced to dynamically adjust the search range. After each iteration, the algorithm checks and constrains the parameter boundaries to ensure the physical rationality and stability of the solution. When the maximum number of iterations is reached or the objective function value converges to a preset threshold, the inversion terminates and the optimal fracture density identification factor result is output. This algorithm ensures that the inversion process can fully explore a broad parameter space to avoid getting trapped in local extrema, while also taking into account the physical constraints of actual seismic data, thus improving the reliability and adaptability of the inversion results. Finally, three-dimensional visualization analysis is performed on the multiple sets of fracture density identification factor data obtained from the inversion. The well profile of the first set of fracture density identification factors is shown below. Figure 11 As shown, the second group of fracture density identification factors is connected to the well profile as follows. Figure 12 As shown, when constructing the rock physics model, the first group of cracks are low-angle cracks, and the second group of cracks are high-angle cracks.

[0255] Figure 11-12 The results show that the profile covers a horizontal distance of 2 km and a time depth of 800-1110 ms. The locations of wells A, B, and C are marked in the figure. Combined with... Figure 6 The fracture density prediction curves from Well B show that in the tight sandstone reservoir section of the quartz development zone at approximately 1330-1362m in the target layer, both sets of fracture densities exhibit a "low-high-low" variation characteristic. Correspondingly, in Figure 11-12 The profile results at well B showed a "low-high-low" anomaly distribution within the same layer, indicating a high degree of vertical fracture development in the area, which is highly consistent with the well inversion results. In contrast, wells A and C did not show significant high-value anomalies in the target layer, indicating weaker fracture development in this area. Overall, the inversion results are in good agreement with the logging data. The fracture density in the target layer of gas well B is significantly higher than that of dry wells A and C, and is consistent with the distribution characteristics of high-permeability areas in the logging interpretation. This result provides crucial information for identifying tight sandstone and potential permeable zones in the study area.

[0256] Figure 13-14These are planar maps of the first and second groups of fracture density identification factors in the target layer of the study area, respectively. The results demonstrate the planar distribution characteristics of the two groups of fracture density identification factors within the study area. In the planar maps, the anomalies at well locations A, B, and C are... Figure 11 and Figure 12 The results of the profile analysis are consistent, further verifying the effectiveness of the proposed method. The high-value areas of the fracture identification factor and the high-permeability areas highly overlap, indicating that the proposed pre-stack seismic inversion method based on multiple sets of fracture density identification factors has good applicability in predicting the fracture distribution in complex fractured tight sandstone gas reservoirs. It can quantitatively evaluate the development intensity and spatial distribution of multiple sets of fractures in tight sandstone gas reservoirs and optimize the sweet spots for natural gas exploration.

[0257] The partially saturated petrophysical model with high-angle multi-scale pore and fracture structure constructed in this invention has excellent adaptability to actual geological reservoirs. The fracture density identification factor obtained by inversion has good consistency with the well logging interpretation results and matches well with the interpreted fracture development section. At the same time, the fracture density identification factor prediction profile obtained by inversion of pre-stack seismic data clearly reveals the spatial distribution of fracture development zone and corresponds well with known gas-producing layers, verifying the effectiveness of the method proposed in this invention in the fine characterization of fractures in tight sandstone gas reservoirs.

[0258] The above description is not intended to limit the present invention, nor is the present invention limited to the examples given above. Any changes, modifications, additions, or substitutions made by those skilled in the art within the scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. A pre-stack inversion method for tight sandstone gas reservoirs using multiple sets of fracture density identification factors, characterized in that, Includes the following steps: Step 1: Data Acquisition and Preprocessing Collect pre-stack seismic gathers, well logging data, mineral composition information, and geological stratification data for the study area; Amplitude-preserving denoising, multiple removal, and dynamic correction preprocessing were performed on pre-stack seismic gathers to obtain P-wave amplitude data at different incident angles; environmental correction and depth matching were performed on well logging data. Step 2: Partially saturated rock physics modeling of multi-scale porous and fractured structures. A partially saturated rock physics model with high-angle multi-scale pore structure was constructed by comprehensively utilizing the Voigt-Reuss-Hill average formula, the improved Xu-White model, the differential equivalent medium theory, the Hudson model, the Schoenberg linear slip model, the improved White-Johnson partially saturated fluid model, and the Brown-Korringa anisotropic Gassmann equation. The improved Xu-white model discretizes the pore space into flat, crack-like pores and equiaxed pores to simulate intergranular pores and dissolution pores in the mineral matrix. The equivalent medium modulus of the porous skeleton of dense sandstone is calculated by comprehensively utilizing the Kuster-Toksöz theory. An improved White-Johnson partially saturated fluid model is constructed, which simplifies partially saturated rocks into two types of periodic strata or spherical patches with different fluid saturation, and calculates the equivalent fluid bulk modulus considering the partially saturated fluid effect. Step 3: Multi-set fracture density rock physics inversion based on well logging data Based on actual well logging data and a multi-scale pore structure partially saturated rock physics model, multiple sets of nonlinear inversions of fracture density were constructed. The inversion process introduced a hybrid optimization algorithm that combines particle swarm optimization algorithm with local gradient search. Step 4: Constructing a crack density identification factor based on Poisson impedance The concept of Poisson impedance is extended to fracture density identification. The elastic parameter variation characteristics of a specific fracture group under different fracture density conditions are analyzed. Through parameter optimization driven by rock physics, an identification factor that is sensitive to fracture density of a specific fracture group but relatively insensitive to matrix lithology and pore fluid is constructed. Step 5: Constructing the AVO approximation equation containing multiple sets of crack density identification factors. The constructed crack density identification factor is directly embedded into the linear approximation of the Fatti high-precision AVO reflection coefficient, and the AVO approximation equation containing multiple sets of crack density identification factors is derived to establish a direct relationship between pre-stack seismic response and crack parameters. Step 6: Pre-stack seismic inversion of tight sandstone gas reservoirs based on multiple sets of fracture density identification factors Based on actual pre-stack seismic records and time-domain pre-stack synthetic seismic records, a least-norm 2 fitting is constructed, and a pre-stack seismic inversion objective function is defined. An improved gray wolf optimization algorithm with Gaussian search threshold constraints is proposed, which introduces adaptive Gaussian perturbation, an elite preservation strategy with threshold constraints, and a hybrid Cauchy-Gaussian mutation mechanism. The improved algorithm is embedded into the pre-stack inversion process to ensure that the inversion process can effectively escape local extrema caused by seismic noise and model errors, quickly and stably converge to the vicinity of the global optimum or a high-quality suboptimal solution, and achieve high-precision fitting. Step 7: Comprehensive Interpretation of Crack Development Characteristics A three-dimensional visualization analysis was performed on the multiple sets of fracture density identification factor data obtained by inversion. Combined with the results of geological structural interpretation, the development degree and spatial distribution pattern of multiple sets of fractures in tight sandstone gas reservoirs were quantitatively evaluated.

2. The pre-stack inversion method for tight sandstone gas reservoirs using multiple sets of fracture density identification factors according to claim 1, characterized in that, Step two involves the physical modeling of partially saturated rock with multi-scale porous structures: 2.1 Construction of the mineral matrix of dense sandstone and calculation of its equivalent modulus The equivalent bulk modulus of the mineral matrix was calculated using the Voigt-Reuss-Hill average. and equivalent shear modulus ,formula: ； ； in, This is the equivalent bulk modulus of the mineral matrix; This is the equivalent shear modulus of the mineral matrix; This represents the volume fraction of clay minerals. This represents the bulk modulus of clay minerals. The shear modulus of clay minerals; This represents the volume fraction of quartz. The bulk modulus of quartz; This is the shear modulus of quartz; 2.2 Construction of Porous Framework and Calculation of Equivalent Modulus in Tight Sandstone The improved Xu-White model discretizes the pore space into flat, fracture-like pores and equiaxed pores to simulate intergranular pores and dissolution pores in the mineral matrix. The equivalent medium modulus of the porous framework in tight sandstone is calculated using the Kuster-Toksöz theory. The formula is: ； ； ； in, This is the equivalent bulk modulus of the mineral matrix; It is the equivalent bulk modulus of the porous skeleton; This is the equivalent shear modulus of the mineral matrix; Porosity refers to the porosity of crack-like pores; Geometric factors related to crack-like pores; Porosity refers to the porosity of equiaxed pores; Geometric factors related to equiaxed pores; This is the equivalent shear modulus of the porous skeleton. The polarization factor associated with crack-like pores; The polarization factor associated with equiaxed pores. These are intermediate variables in the calculation process; 2.3 Construction of the microfractured framework of dense sandstone and calculation of its equivalent modulus The equivalent medium theory of differential media is used to construct a dense sandstone microfractured framework by adding randomly distributed microfractures into a porous framework. The formula for calculating the equivalent medium modulus is as follows: ； ； in, The equivalent bulk modulus of the microcracked skeleton; This represents the volume fraction of microcracks. and A tensor characterizing the influence of microcrack shape; The equivalent bulk modulus of the microcrack; It is the equivalent bulk modulus of the porous skeleton; The equivalent shear modulus of the microcracked skeleton; The equivalent shear modulus of the microcrack; This is the equivalent shear modulus of the porous skeleton. 2.4 Dry rock structure of dense sandstone and calculation of equivalent stiffness matrix For multiple fracture systems, the influence of fractures on anisotropic elasticity is described by combining the Hudson model and the Schoenberg linear slip model. A set of high-angle fractures and a set of low-angle fractures are embedded in a dense sandstone skeleton with a microfracture framework as the background medium to construct a dry dense sandstone rock, in which the background medium is an isotropic medium. For the Group of cracks, their crack density expression: ； in, Crack density; For the first Porosity of the crack group; For the first The average aspect ratio of the group of cracks; The first group of cracks is defined as low-angle cracks, which are rotationally symmetric cracks arranged in a direction perpendicular to the X3 axis. The additional stiffness perturbation caused by the presence of low-angle cracks is given using the Hudson model. The expression for the first-order correction term: ； ； in, The stiffness matrix of dry rock with low-angle cracks; The stiffness matrix of the background medium. The stiffness matrix is ​​the stiffness matrix for the first-order correction term of the added stiffness perturbation; The porosity of the first group of cracks; and Lamé constant of the background medium; , and This is the crack interaction term, used to characterize the crack state; , , , , and The Kronecker function; , , and The unit vector component perpendicular to the crack surface; The second group of cracks is defined as high-angle cracks, which are rotationally symmetric cracks arranged in a direction perpendicular to the X1 axis. They are superimposed using the Schoenberg linear slip model. For the second group of cracks, its normal compliance... and tangential compliance The expression: ； ； in, The crack density is for group 2. The equivalent bulk modulus of the microcracked skeleton; and The correlation coefficient for the crack infill material; The equivalent shear modulus of the microcracked skeleton; The overall flexibility of the fractured dry rock is obtained by adding the flexibility tensor of the background medium to the additional flexibility tensors of multiple sets of fractures. Then, the inverse of this flexibility tensor is used to obtain the equivalent stiffness matrix of the dense sandstone dry rock. ,expression: ； in, Additional compliance for the second group of cracks; 2.5 Characterization of Partially Saturated Fluid Effects By incorporating fracture density parameters into the partially saturated model, an improved White-Johnson partially saturated fluid model is constructed. This model simplifies partially saturated rock into two types of periodic strata or spherical patches with different fluid saturations. The equivalent fluid bulk modulus considering the partially saturated fluid effect is calculated, expressed as: ； ； ； in, Gas saturation; The bulk modulus of the gas; The bulk modulus of water; The equivalent shear modulus of dry dense sandstone; It is a relaxation function; This is the Gassmann low-frequency limiting modulus; The relaxation strength; The imaginary unit; Angular frequency; The modified characteristic relaxation time is used to account for the effect of crack density. The characteristic relaxation time; The coupling coefficient is related to the crack direction; Crack density; 2.6 Modeling and Calculation of Equivalent Stiffness Matrix for Partially Saturated Rock with Multi-Scale Pore and Fracture Structure in Tight Sandstone Using the Brown-Korringa anisotropic Gassmann equations for fluid substitution, the anisotropic equivalent stiffness matrix of a multi-scale porous partially saturated rock model is obtained, expressed as: ； in, The equivalent stiffness matrix of a partially saturated rock model with multi-scale porous and fractured structure in dense sandstone; The equivalent stiffness matrix of dense sandstone dry rock; and The equivalent stiffness matrix components of dense sandstone dry rock; This is the equivalent bulk modulus of the mineral matrix; and The Kronecker function; Porosity; The equivalent fluid bulk modulus that takes into account the effects of partially saturated fluid; This is the equivalent stiffness matrix of the mineral matrix.

3. The pre-stack inversion method for tight sandstone gas reservoirs using multiple sets of fracture density identification factors according to claim 1, characterized in that, Step 3: The process of multi-set fracture density rock physics inversion based on well logging data: Using the model as the forward modeling operator and the P-wave and S-wave velocities from actual well logging data as constraints, the minimum fitting error between the measured well logging data and the model predictions is calculated. Then, multiple sets of fracture densities at each depth point are calculated, thereby inverting multiple sets of fracture density curves along the well depth. The inversion objective function is defined as the error between the measured well logging data and the model predictions. The expression is as follows: in, Forward modeling of the longitudinal wave velocity in rock physics; For logging actual P-wave velocity; The transverse wave velocity is simulated in the forward modeling of the rock physics model; For logging actual shear wave velocity; For regularization terms; The regularization coefficient is used. The crack density is for group 1. The crack density is for group 2. A hybrid optimization algorithm combining particle swarm optimization and local gradient search is used for inversion solution. This algorithm combines the fast convergence characteristics of particle swarm optimization and the ability of local gradient search to escape local optima.

4. The pre-stack inversion method for tight sandstone gas reservoirs using multiple sets of fracture density identification factors according to claim 1, characterized in that, Step four involves the construction of the crack density identification factor based on Poisson impedance: Poisson impedance The general expression: ； ； ； in, For longitudinal wave impedance; For transverse wave impedance; The rotation factor; For longitudinal wave velocity; The transverse wave velocity; For density; based on a rock physics model, the elastic response characteristics of multiple sets of fractures under different fracture densities are analyzed, and the expression of the fracture density sensitivity function is defined: ； in, This refers to a certain elastic parameter under different crack densities in multiple sets of cracks; Crack density; The optimal rotation coefficient is obtained by fitting the maximum correlation coefficient between the crack density sensitivity function and the Poisson impedance. ,expression: ； in, It is a Poisson impedance; It is a crack density sensitive function; For two sets of cracks with different orientations, the concept of Poisson impedance is extended to crack density identification, defining two sets of crack density identification factors, expressed as: ； ； in, The crack density identification factor for the first group of cracks; The optimal rotation coefficients for group 1; The crack density identification factor for the second group of cracks; The optimal rotation coefficients for group 2; For longitudinal wave impedance; This is the transverse wave impedance.

5. The pre-stack inversion method for tight sandstone gas reservoirs using multiple sets of fracture density identification factors according to claim 1, characterized in that, Step five involves the construction of the AVO approximation equation containing multiple sets of crack density identification factors: The two sets of crack density identification factors based on Poisson impedance are extended to a differential expression, as shown in the following expression: ； ； in, This represents the relative change in the crack density identification factor of the first group of cracks at the reflective interface. This represents the relative change in the crack density identification factor of the second group of cracks at the reflective interface. The relative rate of change of the crack density identification factor of the first group of cracks at the reflective interface; The relative rate of change of the crack density identification factor for the second group of cracks at the reflective interface; This represents the relative change in longitudinal wave impedance; This represents the relative change in transverse wave impedance; The optimal rotation coefficients for group 1; The optimal rotation coefficients for group 2; This represents the P-wave velocity; the current value is the same as the P-wave velocity measured in the well logging. same; This represents the shear wave velocity; the current value is the same as the shear wave velocity measured in the well logging. same; Define intermediate variables , , , ,expression: ； ； ； ； The simplified differential expression of the two sets of crack density identification factors based on Poisson impedance is as follows: ； ； By solving the differential expressions for the two sets of crack density identification factors, the relative rate of change of the longitudinal wave impedance is obtained. The relative rate of change of transverse wave impedance ,expression: ； ； Substituting the differential expressions of the two sets of crack density identification factors into the linear approximation of the Fatti high-precision AVO reflection coefficient, an AVO approximation equation containing multiple sets of crack density identification factors is derived. This yields the quantitative relationship between the PP wave reflection coefficient and the crack density identification factor, expressed as: ； in, The angle of incidence is denoted as .

6. The pre-stack inversion method for tight sandstone gas reservoirs using multiple sets of fracture density identification factors according to claim 1, characterized in that, In step six, during the direct inversion of pre-stack earthquakes, an objective function is constructed based on minimizing the fitting error between the actual pre-stack earthquake records and the synthetic earthquake records. The expression is: ； in, This is an actual earthquake record. For theoretically synthesized seismic records, N represents the number of incident angles.

Images