Gas storage geologic body stability intelligent prediction method considering configuration heterogeneity

Abstract

The invention relates to the technical field of reservoir assessment, and discloses a gas storage geologic body stability intelligent prediction method considering configuration heterogeneity, and the method comprises the following steps: 1, reservoir configuration unit division and feature extraction; 2, carrying out coupling analysis on geomechanical parameters and a crustal stress field; step 3, carrying out correlation modeling on configuration heterogeneity and mechanical characteristics; and 4, training and optimizing the stability intelligent prediction model. According to the method, deep coupling characterization of geologic configuration heterogeneity and mechanical response is achieved, and the physical basis and precision of stability prediction are remarkably improved.

Description

Technical Field

[0001] This invention relates to the field of reservoir assessment technology, and in particular to an intelligent prediction method for the stability of gas-bearing geological bodies that takes into account structural heterogeneity. Background Technology

[0002] In the energy sector, including oil and gas exploration and development, reservoir stability prediction is of paramount importance. It involves a comprehensive consideration of reservoir geological characteristics, mechanical properties, and various complex geological conditions. Accurate prediction of reservoir stability provides crucial information for resource extraction planning and ensuring the safety of production operations, directly impacting the efficiency and effectiveness of energy development.

[0003] A search revealed that Chinese patent CN120258045A discloses a method and system for constructing a reservoir evaluation model based on deep learning. This method involves building a specific neural network model, training it with input well logging data, and then predicting features such as reservoir porosity. These methods typically utilize convolutional neural networks to extract features from the input data, and then output the prediction results through traditional neural networks.

[0004] However, existing technologies have significant drawbacks. They suffer from three main limitations: First, they are disconnected from geological configuration and mechanical analysis, often employing homogeneous models or simplified layered models that fail to characterize the differentiated mechanical features of different configuration units such as channels and lobes, resulting in a lack of solid physical foundation for stability predictions and limited accuracy. Second, they lack interpretability and generalization ability; purely data-driven "black box" models lack clear physical rules and constraints, making it difficult to explain the geological mechanisms behind the predictions, and exhibiting poor generalization ability in complex heterogeneous reservoirs. Third, they have low sensitivity to dynamic responses, relying primarily on static geological and mechanical data without effectively integrating dynamic information such as inter-well tracer migration, failing to reflect the coupling relationship between actual reservoir fluid movement and mechanical stability, leading to discrepancies between risk identification and actual reservoir behavior, and making it difficult to meet engineering application requirements.

[0005] To address this, we propose an intelligent prediction method for the stability of gas-bearing geological bodies that takes into account the heterogeneity of their configuration. Summary of the Invention

[0006] The present invention mainly addresses the technical problems existing in the prior art and provides an intelligent prediction method for the stability of gas storage geological bodies that takes into account the heterogeneity of the configuration.

[0007] To achieve the above objectives, the present invention adopts the following technical solution: an intelligent prediction method for the stability of gas-bearing geological bodies considering structural heterogeneity, comprising the following steps:

[0008] Step 1: Reservoir configuration unit division and feature extraction

[0009] Based on core observations, well logging curves, and seismic interpretation data of the target gas storage, the reservoir configuration unit types are identified; the spatial distribution of configuration units is divided, and the geometric, compositional, and structural characteristic parameters of each unit are extracted; a configuration unit classification map and a characteristic parameter dataset are output.

[0010] Step 2: Coupled analysis of geomechanical parameters and geostress field

[0011] For core samples of each configuration unit, rock mechanical parameters were obtained through triaxial compressive strength testing and elastic modulus determination; combined with downhole in-situ stress monitoring data, in-situ stress distribution models for different configuration units were constructed; and a mechanical parameter-in-situ stress coupled dataset was output.

[0012] Step 3: Modeling the correlation between configurational heterogeneity and mechanical characteristics

[0013] Spatially match the characteristic parameter dataset from step 1 with the mechanical parameter-soil stress coupling dataset from step 2; establish a quantitative relationship between configuration characteristic parameters and mechanical parameters through multiple regression analysis.

[0014] Compressive strength = f(clay mineral content, particle sorting coefficient, interlayer density)

[0015] Elastic modulus = g(quartz content, thickness, contact relationship type)

[0016] Identify the associated deviation samples caused by parameter mutations at the boundary of configuration units, correct the quantitative relationship using the weighted least squares method, and generate a configuration-mechanical response association rule library;

[0017] Step 4: Training and Optimization of the Stability Intelligent Prediction Model

[0018] The model uses configuration feature parameters and geostress values ​​as input layers and stability evaluation indicators as output layers. A random forest algorithm is used to train the prediction model, and the model weights are initialized by constraining the configuration-mechanical response association rule base. Dynamic connectivity tracer data is introduced as a regularization term for model training to optimize prediction accuracy. The configuration feature parameters of the target gas storage are input into the trained prediction model, and a stability risk zoning map is output.

[0019] Preferably, step 1, identifying reservoir configuration unit types, specifically includes determining the basic type of configuration unit based on sedimentary structural differences observed in core samples, quantitatively classifying by well logging curves and verifying through seismic interpretation, distinguishing sedimentary structural differences by core bedding types, with channel units exhibiting cross-bedding, lobed main units exhibiting graded bedding, lobed lateral units exhibiting wavy bedding, and outer fan terminal units exhibiting horizontal bedding; the quantitative classification by well logging curves uses a combination of natural gamma curves and resistivity curves, classifying by natural gamma value (GR), with sandy units less than 60 API and argillaceous units greater than 60 API. At 80 API, based on resistivity values ​​RT, sandy units are greater than 50 Ω·m, and clayey units are less than 20 Ω·m. Seismic interpretation and verification were conducted by identifying seismic reflection characteristics: strong strip-shaped reflections indicate waterway units, continuous floret-shaped reflections indicate floret-leaf main body units, weak reflections indicate floret-leaf lateral edge units, and low-amplitude chaotic reflections indicate outer fan terminal units. The final determined configuration unit types include waterways, floret-leaf main bodies, floret-leaf lateral edges, and outer fan terminals. In terms of spatial distribution, each type of unit satisfies the following conditions: waterways are distributed around the floret-leaf main body, floret-leaf lateral edges are located at the edge of the floret-leaf main body, and outer fan terminals are distributed on the outermost periphery.

[0020] Preferably, step 1, which extracts feature parameters, specifically includes geometric feature parameter extraction, compositional feature parameter extraction, and structural feature parameter extraction; among the geometric feature parameters, thickness H...

[0021] The depth difference H = D between the top and bottom interfaces of the sandstone in the well logging curve. 顶 -D 底 D 顶 and D 底 The distribution range S is derived from well logging curve interpretation results; the distribution range S is calculated using the boundary coordinates of the configuration unit from seismic interpretation, where S = ∫∫dxdy, and the integration region is the area enclosed by the boundary of the configuration unit; the contact relationship type T is assigned according to the interface contact characteristics, with scour contact T = 1 and gradual contact T = 0; the quartz content Q in the compositional characteristic parameters is the quartz mass percentage from whole-rock X-ray diffraction analysis, Q = m q / m 总 m q For the mass of quartz, m 总 The total mass of the rock is represented by C; the clay mineral content C is the percentage of clay minerals determined by X-ray diffraction analysis, C = m. c / m 总 m c For clay mineral content, the particle size d in the structural characteristic parameters is the median particle size measured by a laser particle size analyzer, d = d 50 The sorting coefficient σ is calculated using the Fock-Ward formula. φ 84 and φ 16 These are the particle size values ​​corresponding to 84% and 16% in the particle size accumulation curve.

[0022] Preferably, step 2, constructing the geostress distribution model, specifically includes assigning parameters to the finite element model for geostress component calculation, solving the model, and calculating the maximum horizontal principal stress σ in the geostress component calculation. H =k H ×σ v Minimum horizontal principal stress σ h =k h ×σ v k H and k h The horizontal stress coefficient is derived from the inversion results of geostress monitoring, and also includes the vertical stress σ. v σ is calculated using downhole monitoring data. v =ρgh, where ρ is the rock density, g is the gravitational acceleration, and h is the burial depth; the elastic modulus E and Poisson's ratio μ of each configuration element in the finite element model parameter assignment are assigned according to experimental data, specifically: waterway element E = 68.11 GPa, μ = 0.25; main leaf element E = 45.98 GPa, μ = 0.25; leaf lateral edge element E = 37.84 GPa, μ = 0.21; outer fan tip element E = 23.68 GPa, μ = 0.37; the model solution uses finite element software to calculate the geostress distribution M = (σ H ,σ h ,σ v ).

[0023] Preferably, step 3, spatial matching, specifically includes mesh generation, data mapping, and missing value imputation. Mesh generation uses a three-dimensional mesh system with mesh node coordinates (x...). i ,y j ,z k ), i = 1..n, j = 1..m, k = 1..p; Data mapping maps the feature parameters P(x,y,z) of step 1 and the mechanical parameters M(x,y,z) of step 2 to the mesh nodes P(x,y,z) through coordinate matching. i ,y j ,z k )=P1,M(x i ,y j ,z k ) = M1; Missing value imputation uses isomorphic element interpolation. For a missing node (x0, y0, z0), find the k nearest valid nodes within its corresponding isomorphic element and calculate P(x0, y0, z0) = ∑(P i ×w i ) / ∑w i w i The weighting factor is inversely proportional to the distance; the final coupled database DB contains the feature parameters and mechanical parameters of each grid node. DB = {P(x i ,y j ,zk ),M(x i ,y j ,z k )}.

[0024] Preferably, establishing the quantitative relationship in step 3 specifically includes constructing a compressive strength model and an elastic modulus model, with the compressive strength model being σ. c = a×C+b×σ+c×D, where a, b, and c are regression coefficients calculated using multiple linear regression, and σ... c For compressive strength, C is the clay mineral content, σ is the sorting coefficient, and D is the interlayer density; the elastic modulus model is E = d × Q + e × H + f × T, where d, e, and f are regression coefficients, E is the elastic modulus, Q is the quartz content, H is the thickness, and T is the contact relationship type; the regression calculation uses the least squares method with the objective function min∑(measured value - predicted value). 2 The coefficient values ​​are obtained through iterative solving to generate quantitative relationship equations for subsequent construction of the association rule base.

[0025] Preferably, step 3, which corrects the quantitative relationship, specifically includes weighted correction based on deviation sample identification. Deviation sample identification is achieved through residual calculation, where residual e = measured value - predicted value. A deviation sample is marked when |e| > 15% × measured value. Weighted correction employs a weighted least squares method with the objective function min∑(w...) i ×e 2 ), w i Weighting factor w i =1 / (e 2 The biased samples are given greater weight; the corrected compressive strength model is σ. c '=a′×C+b′×σ+c′×D, the elastic modulus model is E′=d′×Q+e′×H+f′×T, a′, b′, c′, d′, e′, f′ are the corrected coefficients obtained by recalculation; all the corrected regression equations are organized into a configuration-mechanical response association rule base, which is output as prior knowledge to step 4.

[0026] Preferably, the model framework construction in step 4 specifically includes an input layer, a hidden layer, and an output layer; the input layer contains configuration feature parameters and geostress value input vector X = [H,σ,D,C,σ]. H ,σ h ,σ vThe hidden layers are set to 3 layers. The first layer has 12 neurons and the activation function is ReLU. The output is h1 = ReLU(W1×X+b1), where W1 is the weight matrix and b1 is the bias term. The second layer has 8 neurons and the output is h2 = ReLU(W2×h1+b2), where W2 and b2 are the parameters of the second layer. The third layer has 4 neurons and the output is h3 = ReLU(W3×h2+b3), where b3 is the parameter of the third layer. The output layer outputs the stability evaluation indicators crack propagation rate v and rock mass deformation u. The output vector is Y = [v,u] = W4×h3+b4, where W4 and b4 are the parameters of the output layer. The model framework initializes the weight matrix W1-W4 through the association rule base so that the initial weights conform to the configuration-mechanical association law.

[0027] Preferably, step 4, training the prediction model, specifically includes sample partitioning model training and regularization optimization; the sample partitioning divides the coupled database DB into a 70% training set and a 30% validation set, with the training set used for model parameter learning and the validation set used for accuracy evaluation; the model training uses a random forest algorithm to construct 100 decision trees, and the splitting features of each tree are selected based on the Gini coefficient: Gini = 1 - ∑(p i ) 2 p i The splitting threshold is determined by the proportion of sample categories, with reference to the association rule base constraint; the regularization optimization introduces dynamic connectivity tracer data to measure the inter-well tracer migration rate v. 示踪 loss function as a regularization term λ is the regularization coefficient, v 预测 The connectivity rate predicted by the model; the coefficient of determination R for the training iteration to the validation set. 2 Stop when ≥0.9 Generate the trained prediction model.

[0028] Preferably, step 4 outputs a stability risk zoning map, specifically including risk level classification, spatial overlay, and visualization output. The risk level classification is based on the predicted crack propagation rate v and rock mass deformation u. Low-risk areas satisfy v < 0.1 mm / d and u < 50 μm; medium-risk areas satisfy 0.1 mm / d ≤ v < 0.5 mm / d and 50 μm ≤ u < 100 μm; high-risk areas satisfy v ≥ 0.5 mm / d and u ≥ 100 μm. Spatial overlay superimposes the risk level and the configuration unit classification map in three-dimensional space. High-risk areas are highlighted in the transition zone between the outer fan tip and the leaf lateral edge, as well as areas with high clay content. The visualization output uses a raster data format, assigning a risk level value R to each grid unit. R = 1, 2, 3 correspond to low, medium, and high risks, generating a stability risk zoning map that includes the spatial distribution of risk levels and statistical values ​​of configuration characteristic parameters for each level, providing an intuitive basis for evaluating the stability of the gas storage facility.

[0029] Beneficial effects

[0030] This invention provides an intelligent prediction method for the stability of gas-bearing geological bodies that considers structural heterogeneity. It has the following beneficial effects:

[0031] (1) This intelligent prediction method for the stability of gas-bearing geological bodies considering configurational heterogeneity achieves deep coupling characterization of geological configurational heterogeneity and mechanical response, significantly improving the physical basis and accuracy of stability prediction. Existing technologies often separate geological description from mechanical analysis or use overly simplified homogeneous models. This scheme achieves deep coupling through three core technologies: First, by integrating core bedding identification, quantitative classification of logging curves (GR / RT threshold), and interpretation of seismic reflection characteristics, a fine classification system is established that includes four types of configurational units and their spatial relationships: waterways, main body of lobes, lateral edge of lobes, and outer fan tip. A dataset of seven structured feature parameters covering geometry, composition, and structure is extracted. Second, when constructing the geostress distribution model, mechanical parameters based on triaxial experiments are assigned according to the differences in configurational unit types. Based on stress component calculation and finite element simulation, a three-dimensional stress tensor that strictly matches the configurational unit is generated. Finally, by using 3D mesh spatial matching and isomorphic element interpolation, the feature parameter set and the mechanical-stress dataset were precisely fused into a coupled database, providing a spatially consistent and attribute-complete data foundation for subsequent modeling. This end-to-end design, from identification and parameterization to spatial coupling, ensures that the characteristics of heterogeneous configurations are directly mapped to mechanical response analysis.

[0032] (2) This intelligent prediction method for the stability of gas-bearing geological bodies considering configurational heterogeneity constructs a hybrid "white box + black box" modeling framework that combines physical interpretability and data adaptability, enhancing the model's generalization ability and reliability. Existing pure data-driven models are often considered "black boxes" lacking physical constraints; while pure mechanistic models struggle to characterize complex heterogeneity. This scheme adopts a two-stage hybrid modeling strategy: In the first stage, explicit configuration-mechanistic quantitative relationships are established through multiple linear regression, and residual-driven error correction is performed using weighted least squares, forming a "white box" association rule base with clear physical meaning. In the second stage, when training the random forest intelligent prediction model, this rule base is used as prior knowledge to constrain model weight initialization, ensuring that the initial response direction of the model conforms to geomechanical laws. At the same time, dynamic connectivity tracer data is introduced as a regularization term, forcing the model to learn actual fluid connectivity characteristics. This design, which uses "white box" rules to guide initialization, "black box" algorithms to learn complex relationships, and additional measured data regularization, improves the ability to handle nonlinear and heterogeneous problems while maintaining physical consistency.

[0033] (3) This intelligent prediction method for the stability of gas-bearing geological bodies considering configurational heterogeneity integrates multi-source dynamic monitoring data to drive model optimization, significantly improving the sensitivity of risk identification to actual reservoir behavior. Existing methods mostly rely on static geological and mechanical data, which are difficult to reflect the dynamic response of the reservoir. This scheme innovatively integrates the dynamic connectivity index of inter-well tracer migration rate into the regularization term of model training. During random forest training, by incorporating the absolute error between the measured tracer migration rate and the connectivity rate predicted by the model into the loss function and optimizing the regularization coefficient, the model must simultaneously match the actual fluid flow characteristics while fitting the relationship between static parameters and target variables. This technique directly establishes a quantitative link between stability risk and reservoir dynamic connectivity, ensuring that the predicted high-risk areas not only reflect static configurational weaknesses but also capture potential seepage-mechanical coupling failure paths revealed by actual fluid movement, greatly improving the engineering practicality of risk warning. Attached Figure Description

[0034] To more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings in the following description are merely exemplary, and those skilled in the art can derive other embodiments based on the provided drawings without creative effort.

[0035] Figure 1 This is a flowchart of the steps of the present invention;

[0036] Figure 2 This is a schematic diagram showing the reservoir configuration unit types and spatial relationships of the present invention;

[0037] Figure 3 This is a flowchart of the coupled analysis of geomechanical parameters and geostress field in this invention;

[0038] Figure 4 This is a schematic diagram of the configuration-mechanical response correlation modeling process of the present invention;

[0039] Figure 5 This is a flowchart illustrating the structure and training process of the intelligent stability prediction model of this invention. Detailed Implementation

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

[0041] Example: An intelligent prediction method for the stability of gas-bearing geological bodies considering configurational heterogeneity, such as... Figures 1-5 As shown, it includes the following steps:

[0042] S1. Reservoir configuration unit division and feature extraction:

[0043] Based on core observations, logging curves (including natural gamma and resistivity curves), and seismic interpretation data of the target gas reservoir, the reservoir configuration unit type is identified.

[0044] In the reservoir configuration unit identification process of step 1, sedimentary structural information is first obtained through core observation, and the basic type of configuration unit is identified based on different bedding characteristics. Specifically, when a large number of cross-beddings appear in the core, it indicates that the area has strong hydrodynamic conditions, indicating a channel unit; if graded bedding is developed, showing a transition from coarse to fine sedimentary grain size, it indicates a leaf-shaped main unit; when wavy bedding is observed, it represents a weakly disturbed sedimentary environment, and is identified as a leaf-shaped lateral unit; if the bedding is mainly parallel horizontal, it usually appears in a distant low-energy environment, corresponding to an outer fan terminal unit. Based on the above core bedding type identification results, the configuration units are further quantitatively divided by well logging curves, with the natural gamma value GR and resistivity value RT as key parameters: GR represents the content of natural radioactive elements in the formation, mainly reflecting the clay content. If GR < 60 API, it indicates that the interval is dominated by sand and belongs to a sand body unit; if GR > 80 API, the clay content is high and it belongs to a clay body unit; intermediate values ​​are determined by resistivity. RT reflects the resistivity of rock strata. Sand bodies typically have higher resistivity due to their higher porosity and lower conductivity of conductive minerals. If RT > 50 Ω·m, the stratum is a sandy unit; if RT < 20 Ω·m, it indicates a clayey unit. This combined judgment method improves the accuracy of structural unit identification. Subsequently, 3D seismic data was used for interpretation and verification. By analyzing reflection characteristics, the spatial consistency of structural identification was enhanced. Channel units exhibit strong amplitude banded reflections in the seismic profile, indicating their good continuity and high density properties. Lobe-shaped main units show wide-amplitude continuous reflections, reflecting their lateral distribution and internal homogeneity. Lobe-shaped lateral edge units show weak reflection bands with strong transitions and unclear boundaries. Outer fan terminal units exhibit low-amplitude, chaotic reflection characteristics, indicating that their deposits have small grain sizes and spatial discontinuities. By combining core observation, quantitative classification of well logging parameters, and seismic reflection characteristics, a complete system of morphological unit types was established, including four categories: channels, main body of lobes, lateral margins of lobes, and outer fan tips. Their spatial relationships were clearly defined: channel units linearly interweave around the main body of lobes, providing sediment supply channels; lateral margins of lobes are located at the edge of the main body, receiving sediment diffusion; and outer fan tips are located at the outermost periphery of the sedimentary system, forming the sedimentary terminus. This identification process provides a precise and clearly categorized data foundation for subsequent extraction of morphological characteristic parameters and input into stability prediction models.

[0045] The spatial distribution of configuration units is divided, and the geometric feature parameters (including thickness, distribution range, and contact relationship), composition feature parameters (including quartz content and clay mineral content), and structural feature parameters (including particle size and sorting coefficient) of each unit are extracted. The configuration unit classification map and feature parameter dataset are output.

[0046] After identifying the structural units, geometric, compositional, and structural feature parameters need to be extracted from each unit to form a complete feature parameter dataset, providing the basic input for subsequent modeling of the correlation between structural configuration and mechanical response. Geometric feature parameters include thickness, distribution range, and contact relationship type. Thickness H represents the vertical scale of the structural unit, and its value is derived from the depth D of the sandstone top interface in the well logging curve interpretation results. 顶 With bottom interface depth D 底 The difference is determined by the formula H = D. 顶 -D 底 This reflects the development intensity of the configuration unit in the vertical direction.

[0047] The extent S measures the spatial extension of a unit in a plane, obtained by integrating the two-dimensional coordinates of the morphological unit boundaries in seismic interpretation results. The formula is S = ∫∫d xdy, where the integration region is the area enclosed by the morphological boundary, representing the degree of occupancy of the unit in the sedimentary body. The contact relationship type T is a qualitative classification index of the interface geometry, assigned based on interface morphology records from core observations. Scour contact represents a strong unconformity and is assigned a value of T = 1, while gradual contact indicates a strong sedimentary transition and is assigned a value of T = 0. This parameter influences the interfacial shear characteristics of the mechanical response. Among the compositional characteristic parameters, the quartz content Q represents the proportion of rigid framework grains, derived from whole-rock X-ray diffraction analysis, and is calculated using the formula Q = m q / m 总 , where m q m is the mass of quartz in the sample. 总 The total mass of the sample is represented by quartz content; a high quartz content typically corresponds to a higher elastic modulus. Clay mineral content (C) measures the proportion of compressible components in a rock; its source is also X-ray diffraction analysis, and the formula is C = m... c / m 总 m c This refers to the mineral content of clay; an increase in its content often leads to a decrease in compressive strength. Structural characteristic parameters include particle size and sorting coefficient. The particle size d represents the median diameter of the particle, measured by a laser particle size analyzer, and is taken as d = d... 50 , representing the average size of particle accumulation in this unit, is related to reservoir tightness and fracture propagation.

[0048] The sorting coefficient σ reflects the dispersion of particle size distribution, and the calculation formula is as follows: Where φ84 and φ 16 The sorting coefficient is used to calculate the logarithmic value of the particle size corresponding to the cumulative percentages of 84% and 16% in the cumulative particle size distribution. The larger the sorting coefficient, the more disordered the particle distribution, which may lead to more severe local stress concentration. All extracted parameters form a structured feature parameter dataset P, including seven indicators: H, S, T, Q, C, d, and σ, which are then passed as input to the configuration and mechanical response correlation modeling module in step 3.

[0049] S2. Coupled analysis of geomechanical parameters and geostress field:

[0050] For core samples of each configuration unit, rock mechanical parameters (including compressive strength, elastic modulus, and Poisson's ratio) are obtained through triaxial compressive strength testing and elastic modulus determination. Combined with downhole in-situ stress monitoring data, in-situ stress distribution models for different configuration units are constructed, and mechanical parameter-in-situ stress coupled datasets are output.

[0051] In step 2, the construction of the geostress distribution model aims to characterize the three-dimensional stress state within different configuration units, serving as the basis for subsequent modeling of the correlation between configuration and mechanical characteristics. First, geostress components are calculated, and the maximum horizontal principal stress σ is obtained using downhole monitoring data. H Minimum horizontal principal stress σ h With vertical stress σ v Three components. Vertical stress σ v Determined based on the weight of the overlying strata, the calculation formula is σ. v =ρgh, where ρ is the rock density in kilograms per cubic meter (derived from density logging data), g is the gravitational acceleration, and h is the burial depth of the measuring point in meters (derived from well depth records). The horizontal stress components are obtained through their proportional relationship with the vertical stress, where the maximum horizontal principal stress σ... H =k H ×σ v Minimum horizontal principal stress σ h =k h ×σ v k H and k hThese are the horizontal stress coefficients, representing the degree of anisotropy of in-situ stress, derived from multi-source in-situ stress monitoring inversion results, including methods such as hydraulic fracturing and microseismic wave velocity inversion. After estimating the in-situ stress components, a finite element model is constructed to simulate the stress distribution within the structural units. The model is divided based on the structural unit boundaries extracted in step 1, and different elastic mechanical parameters are assigned to different units. Specific parameter values ​​are based on triaxial compression test data. The channel unit has the highest rigidity, with an elastic modulus E of 68.11 GPa and a Poisson's ratio μ of 0.25; the main leaf unit is the second most rigid unit, with an E of 45.98 GPa and a μ of 0.25; the leaf lateral edge unit has an elastic modulus of 37.84 GPa and a Poisson's ratio of 0.21, reflecting its transitional structural characteristics; the outer fan terminal unit, due to its high mud content and loose structure, is assigned a value of E = 23.68 GPa and μ = 0.37, representing the lowest mechanical stability.

[0052] After inputting the above parameters into the finite element software, the geostress distribution state of each configuration element is solved based on the static equilibrium equations. The output result is the stress tensor M = (σ H ,σ h ,σ v ), representing the stress level in each principal direction in three-dimensional space within each element. Finally, the stress distribution results of all configuration elements are integrated with the corresponding mechanical parameters into a coupled dataset, which is then passed to step 3 for correlation modeling and analysis of configuration heterogeneity and mechanical response.

[0053] S3. Modeling the correlation between configurational heterogeneity and mechanical properties:

[0054] Spatial matching is performed between the feature parameter dataset of step 1 and the mechanical parameter-soil stress coupling dataset of step 2. Through multivariate regression analysis, a quantitative relationship between configuration feature parameters and mechanical parameters is established, and then a configuration-mechanical response association rule base is generated.

[0055] In step 3, spatial matching is a crucial step in coupling configurational characteristic parameters with geostress mechanical parameters, mainly including three parts: 3D mesh generation, data mapping, and missing value imputation. First, mesh generation is performed to construct a unified 3D computational space. Mesh nodes are generated using regular spacing, and the node coordinates are defined as (x...). i ,y j ,z k ), where i = 1..n, j = 1..m, k = 1..p, representing the grid sequence number along the X, Y, and Z directions. This grid system covers the entire gas storage configuration area, achieving spatial discretization. Subsequently, the configuration feature parameters P(x,y,z) extracted in step 1 and the geostress mechanical parameters M(x,y,z) obtained in step 2 are mapped to the above grid system. Based on the principle of coordinate consistency, at each grid node (x... i ,yj ,z k Assign values ​​to P(x) respectively. i ,y j ,z k )=P1,M(x i ,y j ,z k = M1, where P1 and M1 are attribute values ​​sampled or calculated at corresponding spatial locations. P includes thickness H, distribution range S, contact relationship T, quartz content Q, clay mineral content C, particle size d, and sorting coefficient σ. M includes the maximum horizontal principal stress σ. H Minimum horizontal principal stress σ h Vertical stress σ v And elastic modulus E, Poisson's ratio μ.

[0056] For missing values ​​in some nodes due to insufficient data coverage, isomorphic unit interpolation is used to fill them in. Within the isomorphic unit to which the missing node (x0, y0, z0) belongs, the k nearest valid neighboring nodes (x, y0, z0) are found. i ,y j ,z k ), and calculate the feature parameter value of the node as P(x0,y0,z0)=∑(P i ×w i ) / ∑w i , where P i The parameter values ​​for neighboring nodes, w i The weighting factor is inversely proportional to the distance, meaning the closer the distance, the greater the weight, ensuring spatial continuity of the interpolation results in terms of configuration features. After completing the mapping and interpolation of all node data, the final coupled database DB is generated, with the structure DB = P(x i ,y j ,z k ),M(x i ,y j ,z k ), representing the complete set of parameters on each spatial grid node, this database will be passed as input to the next configuration-mechanical response correlation modeling module to build the training dataset for the stability prediction model.

[0057] In step 3, the correlation modeling between configurational heterogeneity and mechanical characteristics is achieved by establishing a quantitative relationship between compressive strength and elastic modulus. This process is based on the configurational characteristic parameters extracted in step 1 and the mechanical parameters obtained in step 2, and two regression models are constructed using the multiple linear regression method. The expression for the compressive strength model is σ. c = a×C+b×σ+c×D, where σ cThe compressive strength of the structural unit is an important indicator for evaluating stability. C is the clay mineral content, which reflects the compressibility and weak cementation of the rock. σ is the sorting coefficient, which represents the uniformity of the particle composition and is sensitive to local stress concentration. D is the interlayer density, which represents the development strength of non-reservoir interlayers and is closely related to the discontinuity and failure surface distribution within the structure. The coefficients a, b, and c are regression weights, which represent the degree of influence of each variable on the compressive strength. The expression for the elastic modulus model is E = d × Q + e × H + f × T, where E is the elastic modulus, reflecting the stress response capability of the configuration unit in the elastic stage; Q is the quartz content, representing the proportion of the rigid skeleton, which is positively correlated with the elastic characteristics of the rock; H is the thickness, describing the vertical scale of the configuration unit; the greater the thickness, the better it can buffer stress disturbances; and T is the contact relationship type. When it is scour contact, T = 1, indicating strong interface rigidity and strong shear force transmission capability; when it is gradual contact, T = 0, indicating strong interface transition and easy formation of weak surfaces. The coefficients d, e, and f correspond to the weights of the three variables, respectively.

[0058] The regression calculation uses the least squares method, by constructing the objective function min∑(measured value - predicted value). 2 The prediction errors of all sample points are summed, and the parameter coefficients are iteratively updated to minimize the objective function, thus obtaining the optimal regression equation. This process is trained on data pairs of each grid node in the spatially matched database to ensure that the model has physical consistency and statistical significance in the relationship between configuration differences and mechanical response. The two sets of quantitative relationship equations are then integrated into the configuration-mechanical response association rule base, serving as constraints and knowledge base inputs to the training and optimization phase of the stability intelligent prediction model, thereby improving the model's ability to discriminate the stability of different configuration units.

[0059] In step 3, after establishing the quantitative relationship, further error correction is needed to improve the model's prediction accuracy. This process includes two operations: biased sample identification and weighted correction. Biased sample identification is accomplished through residual calculation, where the residual is defined as e = measured value - predicted value, representing the degree of difference between the model's prediction and the actual observation. When the absolute value of the residual exceeds 15% of the measured value, i.e., |e|>0.15×measured value, the sample is marked as a biased sample, indicating that the sample was not effectively modeled in the initial regression. This usually occurs in regions with abnormal configuration characteristics and prominent nonlinear responses. To improve the model's ability to fit such samples, weighted least squares regression correction is used, with the objective function being min∑(w i ×e 2 ), where w i Assign a value of w to the sample weighting factor. i =1 / (e 2 The larger the residual, the higher the weight, prompting the model to pay more attention to the fitting quality of biased samples when re-regressing.

[0060] Based on this, the compressive strength model and the elastic modulus model are reconstructed, and the revised compressive strength relationship is σ. c ′=a′×C+b′×σ+c′×D, where σ c ' represents the corrected predicted compressive strength value, and a', b', and c' are the coefficients obtained after weighted regression; the corrected elastic modulus relationship is E'=d'×Q+e′×H+f'×T, where E' is the corrected predicted elastic modulus value, and d′, e', and f' are the coefficients calculated in the new round of regression.

[0061] All correction operations are based on the spatial mapping relationship between the aforementioned configurational characteristic parameters C, σ, D, Q, H, and T and the target mechanical parameters, ensuring the consistency of the model structure. After weighted correction, the root mean square error of the model residuals is significantly reduced to below 8%, indicating that the corrected model has better generalization ability and local adaptability. Finally, all corrected regression equations are compiled into a configuration-mechanical response association rule base, which is output as prior knowledge to step 4 to guide the parameter initialization and weight constraints of the intelligent prediction model, improving the accuracy and interpretability of stability risk identification.

[0062] S4. Training and Optimization of Stability Intelligent Prediction Model:

[0063] The input layer uses structural feature parameters (thickness, sorting coefficient, interlayer density, clay mineral content) and geostress values, while the output layer uses stability evaluation indicators (fracture propagation rate, rock mass deformation). A random forest algorithm is used to train the prediction model, and the model weights are initialized using a configuration-mechanical response association rule base. Dynamic connectivity tracer data (inter-well tracer migration rate) is introduced as a regularization term for model training to optimize prediction accuracy. The structural feature parameters of the target gas storage are input into the trained prediction model, and a stability risk zoning map is output.

[0064] In step 4, the intelligent stability prediction model is constructed using the aforementioned configuration features and geostress parameters as input. A multi-layer neural network structure is used to jointly predict crack propagation rate and rock mass deformation. The model structure includes an input layer, a hidden layer, and an output layer. The input layer receives seven configuration and geostress parameters, forming the input vector X = [H, σ, D, C, σ]. H ,σ h ,σ v ], where H is the thickness of the morphological unit, representing its vertical development intensity; σ is the sorting coefficient, reflecting the uniformity of the particle structure; D is the interlayer density, representing the degree of disturbance to stability by the proportion of non-reservoir layers; C is the clay mineral content, measuring the compressibility and weak cementation of the rock mass; σ H σ h With σ vThese are the maximum horizontal principal stress, minimum horizontal principal stress, and vertical stress, respectively. They are the core indicators for controlling the crack initiation direction and rock mass deformation path. The data are all from the output results of steps 1 and 2 above.

[0065] The model has three hidden layers. The first hidden layer contains 12 neurons, and its output is h1 = ReLU(W1×X+b1), where W1 is the weight matrix from the input layer to the first hidden layer, b1 is the bias term, and the activation function uses ReLU to achieve non-linear transformation.

[0066] The second hidden layer contains 8 neurons, and the output is h2 = ReLU(W2×h1+b2). ​​The third hidden layer contains 4 neurons, and the output is h3 = ReLU(W3×h2+b3). The weight matrices and bias terms of each layer are W2, b2 and W3, b3, respectively. The model parameters are updated through backpropagation optimization.

[0067] The output layer is used to predict stability evaluation indicators. The output vector is Y = [v, u] = W4 × h3 + b4, where v is the crack propagation rate, measuring the dynamic response capability of local rock mass instability, u is the rock mass deformation, representing the overall deformation degree of the configuration unit under in-situ stress, and W4 and b4 are the output layer weights and biases. To ensure that the model has a reasonable physical response trend in the initial stage, the model framework incorporates the results of the configuration-mechanical response association rule library from the previous stage as prior knowledge during construction, initializing the weight matrices W1 to W4 so that they can reflect the basic directionality and sensitivity of the mechanical parameter response when the configuration parameters change, further improving the model's learning efficiency and prediction accuracy in subsequent training.

[0068] In step 4, the training of the stability intelligent prediction model is based on the configuration and geostress coupling database (DB), and completes three core processes sequentially: sample partitioning, model training, and regularization optimization. In the sample partitioning stage, all spatial grid node samples in the coupling database DB are divided into a training set and a validation set, with a ratio of 70% and 30%, respectively. The training set is used to learn the mapping relationship between input parameters and target variables, while the validation set is used to independently evaluate the model's generalization ability. Model training uses the random forest algorithm to construct an ensemble regression model. This algorithm contains 100 decision trees, each of which constructs an independent subset by bootstrapping the training samples, enhancing the model's robustness and nonlinear expressive power. During the feature splitting process at each node, the Gini coefficient is used to select the optimal splitting attribute. The formula for calculating the Gini coefficient is Gini = 1 - ∑(p i ) 2 , where p i This represents the proportion of samples of type i in the current node. This indicator measures the purity of the node. The smaller the Gini value, the more concentrated the samples are and the better the splitting effect.

[0069] Meanwhile, to improve the physical interpretability of feature splitting, the configuration-mechanical response association rule base of step 3 is introduced as a priori constraint to guide the splitting thresholds of key features, such as C, σ, D and σ c The correspondence between these factors can serve as an important basis for determining the trend of fracture propagation rate v. To prevent model overfitting and enhance its sensitivity to actual connectivity characteristics, dynamic connectivity tracer data is introduced as a regularization term. In each training round, the measured inter-well tracer migration rate v is included. 示踪 The connectivity rate v predicted by the model 预测 The absolute error between the two values ​​is incorporated into the loss function, and the overall loss function is defined as follows: Where λ is the regularization coefficient, used to control the weight of the regularization term on the total loss. This term ensures that the model output can reflect the actual changes in reservoir fluid connectivity and improve the ability to identify stability risks under heterogeneous structures.

[0070] The training process continues iteratively until the determination coefficients on the validation set satisfy R. 2 The formula for the coefficient of determination is as follows: up to ≥0.9.

[0071]

[0072] Where ytest represents the true value of the validation set samples, ytest 预测 These are the model's predicted values. Rm is the mean of the measured values. This index reflects the degree to which the predicted values ​​fit the sample fluctuations. 2 The closer the value is to 1, the higher the model accuracy. Ultimately, the model achieves the target R... 2 After the threshold is reached, the trained stability prediction model is output, which is used to generate the stability risk partition map at the end of the step.

[0073] In step 4, the output process of the stability risk zoning map is based on the model prediction results, and sequentially completes three steps: risk level classification, three-dimensional spatial overlay, and visualization generation. First, the stability risk is quantitatively classified according to the crack propagation rate v and rock mass deformation u output by the prediction model. Based on engineering stability control standards, three risk level classification thresholds are set:

[0074] When v < 0.1 mm / d and u < 50 μm, it indicates that the structural unit is dense and the stress disturbance response is low, and it is classified as a low-risk zone; when 0.1 mm / d ≤ v < 0.5 mm / d and 50 μm ≤ u < 100 μm, it indicates that the rock mass has a certain degree of structural heterogeneity and stress sensitivity, and it is classified as a medium-risk zone; when v ≥ 0.5 mm / d and u ≥ 100 μm, it indicates that the crack activity is strong, the deformation is severe, the structural heterogeneity and mechanical weakening are significant, and it is classified as a high-risk zone.

[0075] Subsequently, the risk level of each spatial grid node is superimposed with the configuration unit classification map output in step 1 in three-dimensional space. Spatial coordinate matching is used to achieve a one-to-one spatial correspondence between configuration type and stability level, and the distribution characteristics of high-risk areas in different configuration units are identified.

[0076] Analysis results show that high-risk areas (risk level 3) are mainly concentrated in the outer fan terminal unit and the transition zone between the leaf lateral edge unit and the main unit. These areas generally have high clay mineral content (C), high interlayer density (D), and high sorting coefficient (σ), which leads to low elastic modulus (E) and high compressive strength (σ). c It is lower and more sensitive to ground stress disturbances.

[0077] When finally completing the visualization output, the partitioning results are organized using a raster data format. Each cell in the 3D mesh system is assigned a risk level value R, where R=1 represents low risk, R=2 represents medium risk, and R=3 represents high risk, forming a stability risk partitioning map. The stability risk partitioning map not only shows the distribution of each risk level in 3D space but also simultaneously outputs the statistical results of the main configurational characteristic parameters within each level block, including H, C, D, σ, and σ². H The mean and standard deviation of variables provide a visual reference and quantitative support for subsequent stability control strategies and injection-production scheme optimization, thereby improving the accuracy and foresight of geological body risk management.

[0078] The working principle of this invention is as follows: This scheme achieves deep coupling through three core technologies: First, by integrating core bedding identification, quantitative classification of well logging curves (GR / RT threshold), and interpretation of seismic reflection characteristics, a refined classification system is established, encompassing four types of structural units—channels, main body of lobes, lateral edges of lobes, and outer fan tips—and their spatial relationships. A dataset of seven structured feature parameters covering geometry, composition, and structure is extracted. Second, when constructing the geostress distribution model, mechanical parameters based on triaxial experiments are assigned according to the differences in structural unit types. Based on stress component calculations and finite element simulations, a three-dimensional stress tensor that strictly matches the structural units is generated. Finally, through three-dimensional mesh spatial matching and isomorphic unit interpolation, the feature parameter set P and the mechanical-stress dataset M are precisely fused into a coupled database, providing a spatially consistent and attribute-complete data foundation for subsequent modeling. This end-to-end design, from identification and parameterization to spatial coupling, ensures that heterogeneous structural features are directly mapped to mechanical response analysis.

[0079] This scheme employs a two-stage hybrid modeling strategy: In the first stage, explicit configuration-mechanical quantitative relationships are established through multiple linear regression, and weighted least squares is used for residual-driven error correction, forming a "white-box" association rule base with clear physical meaning. In the second stage, when training the random forest intelligent prediction model, this rule base is used as prior knowledge to constrain model weights during initialization, ensuring that the model's initial response direction conforms to geomechanical laws. Simultaneously, dynamic connectivity tracer data is introduced as a regularization term, forcing the model to learn actual fluid connectivity characteristics. This design—"white-box" rule-guided initialization, "black-box" algorithm learning complex relationships, and additional measured data regularization—improves the ability to handle nonlinear and heterogeneous problems while maintaining physical consistency.

[0080] This innovative approach integrates the dynamic connectivity index of inter-well tracer migration rate into the regularization term for model training. During random forest training, the absolute error between the measured tracer migration rate and the model-predicted connectivity rate is incorporated into the loss function, and the regularization coefficient is optimized. This ensures that the model, while fitting the relationship between static parameters and the target variable, simultaneously matches the actual fluid flow characteristics. This technique directly establishes a quantitative link between stability risk and reservoir dynamic connectivity, ensuring that predicted high-risk areas not only reflect static configuration weaknesses but also capture potential seepage-mechanical coupling failure paths revealed by actual fluid movement, significantly improving the engineering practicality of risk warning.

[0081] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely illustrative of the principles of the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed. The scope of protection of this invention is defined by the appended claims and their equivalents.

Claims

1. A smart prediction method for the stability of gas-bearing geological bodies considering structural heterogeneity, characterized in that, Includes the following steps: Step 1: Reservoir configuration unit division and feature extraction Based on core observations, well logging curves, and seismic interpretation data of the target gas storage, the reservoir configuration unit types are identified; the spatial distribution of configuration units is divided, and the geometric, compositional, and structural characteristic parameters of each unit are extracted; a configuration unit classification map and a characteristic parameter dataset are output. Step 2: Coupled analysis of geomechanical parameters and geostress field For core samples of each configuration unit, rock mechanical parameters were obtained through triaxial compressive strength testing and elastic modulus determination; combined with downhole in-situ stress monitoring data, in-situ stress distribution models for different configuration units were constructed; and a mechanical parameter-in-situ stress coupled dataset was output. Step 3: Modeling the correlation between configurational heterogeneity and mechanical characteristics Spatially match the feature parameter dataset from step 1 with the mechanical parameter-soil stress coupling dataset from step 2; establish a quantitative relationship between configuration feature parameters and mechanical parameters through multiple regression analysis: compressive strength = f(clay mineral content, particle sorting coefficient, interlayer density). Elastic modulus = g(quartz content, thickness, contact relationship type) Identify the associated deviation samples caused by parameter mutations at the boundary of configuration units, correct the quantitative relationship using the weighted least squares method, and generate a configuration-mechanical response association rule library; Step 4: Training and Optimization of the Stability Intelligent Prediction Model The model uses configuration feature parameters and geostress values ​​as input layers and stability evaluation indicators as output layers. A random forest algorithm is used to train the prediction model, and the model weights are initialized by constraining the configuration-mechanical response association rule base. Dynamic connectivity tracer data is introduced as a regularization term for model training to optimize prediction accuracy. The configuration feature parameters of the target gas storage are input into the trained prediction model, and a stability risk zoning map is output.

2. The intelligent prediction method for the stability of gas-bearing geological bodies considering configurational heterogeneity as described in claim 1, characterized in that: Step 1, identifying reservoir configuration unit types, specifically includes determining the basic type of configuration unit based on sedimentary structural differences observed in core samples, quantitatively classifying them using well logging curves and verifying the results through seismic interpretation, distinguishing sedimentary structural differences by core bedding types, identifying channel units as having cross-bedding, lobed main units as having graded bedding, lobed lateral units as having wavy bedding, and outer fan terminal units as having horizontal bedding; quantitative classification using well logging curves employs a combination of natural gamma ray curves and resistivity curves, classifying them according to the natural gamma ray value (GR): sandy units are less than 60 API, and argillaceous units are greater than 80 API. According to the resistivity value RT, sandy units are greater than 50 Ω·m, and clayey units are less than 20 Ω·m. Seismic interpretation and verification are verified by identifying seismic reflection characteristics: strong strip-shaped reflections are waterway units, continuous floret-shaped reflections are floret-leaf main body units, weak reflections are floret-leaf lateral edge units, and low-amplitude chaotic reflections are outer fan terminal units. The final determined configuration unit types include waterways, floret-leaf main bodies, floret-leaf lateral edges, and outer fan terminals. In terms of spatial distribution, each type of unit satisfies the following conditions: waterways are distributed around the floret-leaf main body, floret-leaf lateral edges are located at the edge of the floret-leaf main body, and outer fan terminals are distributed on the outermost periphery.

3. The intelligent prediction method for the stability of gas-bearing geological bodies considering configurational heterogeneity as described in claim 1, characterized in that: The extraction of feature parameters in step 1 specifically includes geometric feature parameter extraction, composition feature parameter extraction, and structural feature parameter extraction; In the geometric characteristic parameters, the thickness H is the difference between the depth of the top and bottom interfaces of the sandstone in the well logging curve, H = D. 顶 -D 底 D 顶 and D 底 The distribution range S is derived from well logging curve interpretation results; the distribution range S is calculated using the boundary coordinates of the configuration unit from seismic interpretation, where S = ∫∫dxdy, and the integration region is the area enclosed by the boundary of the configuration unit; the contact relationship type T is assigned according to the interface contact characteristics, with scour contact T = 1 and gradual contact T = 0; the quartz content Q in the compositional characteristic parameters is the quartz mass percentage from whole-rock X-ray diffraction analysis, Q = m q / m 总 m q For quartz quality, v 总 The total mass of the rock is represented by C; the clay mineral content C is the percentage of clay minerals determined by X-ray diffraction analysis, C = m. c / m 总 m c For clay mineral content, the particle size d in the structural characteristic parameters is the median particle size measured by a laser particle size analyzer, d = d 50 The sorting coefficient σ is calculated using the Fock-Ward formula. φ 84 and φ 16 These are the particle size values ​​corresponding to 84% and 16% in the particle size accumulation curve.

4. The intelligent prediction method for the stability of gas-bearing geological bodies considering configurational heterogeneity as described in claim 1, characterized in that: Step 2, constructing the geostress distribution model, specifically includes calculating the geostress components, assigning parameters to the finite element model, solving the model, and calculating the maximum horizontal principal stress σ in the geostress component calculation. H =k H ×σ v Minimum horizontal principal stress σ h =k h ×σ v k H and k h The horizontal stress coefficient is derived from the inversion results of geostress monitoring, and also includes the vertical stress σ. v σ is calculated using downhole monitoring data. v =ρgh, where ρ is the rock density, g is the gravitational acceleration, and h is the burial depth; the elastic modulus E and Poisson's ratio μ of each configuration element in the finite element model parameter assignment are assigned according to experimental data, specifically: waterway element E = 68.11 GPa, μ = 0.25; main leaf element E = 45.98 GPa, μ = 0.25; leaf lateral edge element E = 37.84 GPa, μ = 0.21; outer fan tip element E = 23.68 GPa, μ = 0.37; the model solution uses finite element software to calculate the geostress distribution M = (σ H ,σ h ,σ v ).

5. The intelligent prediction method for the stability of gas-bearing geological bodies considering configurational heterogeneity as described in claim 1, characterized in that: Step 3, spatial matching, specifically includes mesh generation, data mapping, and missing value imputation. Mesh generation uses a three-dimensional mesh system with mesh node coordinates (x...). i ,y j ,z k ), i = 1..n, j = 1..m, k = 1..p; Data mapping maps the feature parameters P(x,y,z) of step 1 and the mechanical parameters M(x,y,z) of step 2 to the mesh nodes P(x,y,z) through coordinate matching. i ,y j ,z k )=P1,M(x i ,y j ,z k ) = M1; Missing value imputation uses isomorphic element interpolation. For a missing node (x0, y0, z0), find the k nearest valid nodes within its corresponding isomorphic element and calculate P(x0, y0, z0) = ∑(P i ×w i ) / ∑w i w i The weighting factor is inversely proportional to the distance; the final coupled database DB contains the feature parameters and mechanical parameters of each grid node. DB = {P(x i ,y j ,z k ),M(x i ,y j ,z k )}.

6. The intelligent prediction method for the stability of gas-bearing geological bodies considering configurational heterogeneity as described in claim 1, characterized in that: Step 3, establishing the quantitative relationship, specifically includes constructing a compressive strength model and an elastic modulus model. The compressive strength model is σ. c = a×C+b×σ+c×D, where a, b, and c are regression coefficients calculated using multiple linear regression, and σ... c For compressive strength, C is clay mineral content, σ is sorting coefficient, and D is interlayer density; the elastic modulus model is E=d×Q+e×H+f×T, where d, e, and f are regression coefficients, E is elastic modulus, Q is quartz content, H is thickness, and T is contact relationship type. The regression calculation uses the least squares method with the objective function being min∑(measured value - predicted value). 2 The coefficient values ​​are obtained through iterative solving to generate quantitative relationship equations for subsequent construction of the association rule base.

7. The intelligent prediction method for the stability of gas-bearing geological bodies considering configurational heterogeneity as described in claim 1, characterized in that: Step 3, specifically the correction of the quantitative relationship, includes weighted correction based on deviation sample identification. Deviation sample identification is achieved through residual calculation, where residual e = measured value - predicted value. Samples with |e| > 15% × measured value are marked as deviation samples. Weighted correction employs weighted least squares with the objective function min∑(w...). i ×e 2 ), w i Weighting factor w i =1 / (e 2 The biased samples are given greater weight; the corrected compressive strength model is σ. c E' = a'×C+b'×σ+c'×D, the elastic modulus model is E' = d'×Q+e'×H+f'×T, a', b', c', d', e', f' are the corrected coefficients obtained by recalculation; all the corrected regression equations are organized into a configuration-mechanical response association rule base, which is output as prior knowledge to step 4.

8. The intelligent prediction method for the stability of gas-bearing geological bodies considering configurational heterogeneity according to claim 1, characterized in that: Step 4, in constructing the model framework, specifically includes an input layer, a hidden layer, and an output layer; the input layer contains configuration feature parameters and geostress values ​​input vector X = [H, σ, D, C, σ]. H ,σ h ,σ v The hidden layers are set to 3 layers. The first layer has 12 neurons and the activation function is ReLU. The output is h1 = ReLU(W1×X+b1), where W1 is the weight matrix and b1 is the bias term. The second layer has 8 neurons and the output is h2 = ReLU(W2×h1+b2), where W2 and b2 are the parameters of the second layer. The third layer has 4 neurons and the output is h3 = ReLU(W3×h2+b3), where b3 is the parameter of the third layer. The output layer outputs the stability evaluation indicators crack propagation rate v and rock mass deformation u. The output vector is Y = [v,u] = W4×h3+b4, where W4 and b4 are the parameters of the output layer. The model framework initializes the weight matrix W1-W4 through the association rule base so that the initial weights conform to the configuration-mechanical association law.

9. The intelligent prediction method for the stability of gas-bearing geological bodies considering configurational heterogeneity according to claim 1, characterized in that: Step 4, training the prediction model, specifically includes sample partitioning, model training, and regularization optimization. The sample partitioning divides the coupled database DB into a 70% training set and a 30% validation set. The training set is used for model parameter learning, and the validation set is used for accuracy evaluation. The model training employs a random forest algorithm to construct 100 decision trees. The splitting features of each tree are selected based on the Gini coefficient: Gini = 1 - Σ(p...). i ) 2 p i The splitting threshold is determined by the proportion of sample categories, with reference to the association rule base constraint; the regularization optimization introduces dynamic connectivity tracer data to measure the inter-well tracer migration rate v. 示踪 loss function as a regularization term λ is the regularization coefficient, v 预测 The connectivity rate predicted by the model; the coefficient of determination R for the training iteration to the validation set. 2 Stop when ≥0.9 Generate the trained prediction model.

10. The intelligent prediction method for the stability of gas-bearing geological bodies considering configurational heterogeneity according to claim 1, characterized in that: Step 4 outputs a stability risk zoning map, specifically including risk level classification, spatial overlay, and visualization output. The risk level classification is based on the predicted crack propagation rate v and rock mass deformation u. Low-risk areas satisfy v < 0.1 mm / d and u < 50 μm; medium-risk areas satisfy 0.1 mm / d ≤ v < 0.5 mm / d and 50 μm ≤ μ < 100 μm; high-risk areas satisfy v ≥ 0.5 mm / d and u ≥ 100 μm. Spatial overlay superimposes the risk level and configuration unit classification map in three-dimensional space. High-risk areas are highlighted in the transition zone between the outer fan tip and the leaf lateral edge, as well as areas with high clay content. The visualization output uses a raster data format, assigning a risk level value R to each grid cell. R = 1, 2, 3 correspond to low, medium, and high risks, generating a stability risk zoning map that includes the spatial distribution of risk levels and statistical values ​​of configuration characteristic parameters for each level, providing an intuitive basis for evaluating the stability of the gas storage facility.

Images