A method for judging the cracking of supercritical carbon dioxide heating low maturity shale

Abstract

This invention provides a method for determining the initiation of fractures in low-maturity shale heated by supercritical carbon dioxide, belonging to the technical field of low-maturity shale. This invention obtains multi-field coupling parameters by integrating conventional well logging, core experiments, and formation testing data. It then calibrates key coefficients through high-temperature mechanical experiments, hydrocarbon generation pressurization experiments, pyrolysis kinetic experiments, and supercritical carbon dioxide immersion experiments. These parameters are substituted into the multi-field coupling fracture initiation pressure formula to calculate the fracture initiation pressure at each moment. The input to the physical constraint graph attention recursive network outputs corrected fracture initiation pressure and safety margin. When the safety margin is lower than the warning threshold, a two-step approximation method is used to predict the fracture initiation time and adjust the construction parameters. This solves the technical problem of being unable to accurately determine the fracture initiation time of low-maturity shale heated by supercritical carbon dioxide under multi-field coupling conditions in real time.

Description

Technical Field

[0001] This invention belongs to the field of low-maturity shale technology, specifically, it relates to a method for judging the initiation of cracks in low-maturity shale heated by supercritical carbon dioxide. Background Technology

[0002] In-situ refining of medium- and low-maturity shale is an emerging technology for oil and gas resource extraction, which involves injecting supercritical carbon dioxide into the formation and heating it to pyrolyze kerogen and generate hydrocarbons. During engineering implementation, fracture initiation assessment is a crucial step in ensuring construction safety and efficiency. Current technologies typically rely on a single mechanical criterion or empirical formula for predicting fracture initiation pressure, such as the Hubbert-Willis formula based on minimum principal stress or the Kirsch solution based on elasticity. Some approaches incorporate finite element numerical simulations to analyze the temperature or pressure fields separately, while a few studies attempt to incorporate thermal effects into geostress correction. These methods have been widely applied in conventional hydraulic fracturing engineering.

[0003] However, the fracturing process of low-maturity shale heated by supercritical carbon dioxide involves strong coupling between thermal, flow, solid mechanical, and chemical reaction fields. Specifically, this manifests as follows: kerogen pyrolysis generates hydrocarbons, leading to a dynamic increase in pore pressure; supercritical carbon dioxide dissolves carbonate minerals and adsorbs them onto the organic matter surface, resulting in a continuous weakening of rock mechanical strength; and rising formation temperature causes thermal expansion stress superposition. Multiple physicochemical processes are coupled temporally and spatially. Existing single-field or dual-field coupled models cannot simultaneously describe the dynamic evolution of these multiple fields. Empirical formulas lack the ability to quantitatively characterize pyrolysis kinetics and mineral weakening mechanisms, and purely data-driven models have severely insufficient generalization ability under sparse experimental sample conditions, leading to significant deviations between predicted fracturing pressures and actual values.

[0004] In other words, existing technologies have the technical problem of being unable to accurately determine the moment of fracturing in low-maturity shale under supercritical carbon dioxide heating in real time under multi-field coupling conditions. Summary of the Invention

[0005] In view of this, the present invention provides a method for judging the initiation of cracks in low-maturity shale heated by supercritical carbon dioxide, which can solve the technical problem in the prior art that it is impossible to judge the initiation time of low-maturity shale heated by supercritical carbon dioxide in real time and accurately under multi-field coupling conditions.

[0006] This invention is implemented as follows: This invention provides a method for determining the initiation of cracks in low-maturity shale heated by supercritical carbon dioxide, comprising the following steps:

[0007] By integrating conventional well logging, core experiments and formation testing data, the minimum horizontal principal stress, maximum horizontal principal stress, vertical stress, original formation pore pressure, formation temperature, total organic carbon content, carbonate mineral volume fraction, initial tensile strength, elastic modulus, Poisson's ratio, coefficient of thermal expansion, rock density, as well as supercritical carbon dioxide injection pressure, target heating temperature, heating time and heating rate parameters of medium and low maturity shale formations were obtained.

[0008] Based on shale core samples of medium and low maturity, high temperature rock mechanics experiments, high temperature and high pressure closed hydrocarbon generation pressurization experiments, pyrolysis kinetics experiments, and mechanical comparison experiments after supercritical carbon dioxide immersion were carried out. The tensile strength temperature decay coefficient, pressurization conversion coefficient, discrete activation energy distribution model reaction kinetic parameters, strength weakening amplitude coefficient, and strength weakening rate coefficient were calibrated, and the effective porosity and pore throat connectivity coefficient were measured.

[0009] The obtained parameters and calibration parameters are substituted into the multi-field coupled fracturing pressure formula to calculate the fracturing pressure at each time. Then, the physical constraint graph attention recursive network is input to output the corrected fracturing pressure and safety margin. The multi-field coupled fracturing pressure formula integrates the pressurization of kerogen pyrolysis hydrocarbon generation, the intensity reduction of supercritical carbon dioxide, the thermal expansion stress correction and the original geostress into the same analytical framework.

[0010] Substituting the safety margin, the rate of change of the safety margin, and the uncertainty of the initiation pressure into the Ritz adaptive approximation algorithm, the spatial distribution of the initiation pressure and the confidence interval of the initiation pressure at each time point are output.

[0011] Substitute the normalized values ​​of safety margin, safety margin change rate, and fracturing pressure uncertainty into the dynamic adjustment exponential function to calculate the dynamic adjustment exponential. Adjust the time step and iteration strategy of the physical constraint graph attention recursive network according to the interval to which the dynamic adjustment exponential belongs, and update the safety margin according to the difference between the bottom hole pressure and the fracturing pressure obtained from real-time monitoring.

[0012] If the safety margin is lower than the warning threshold, the expected crack initiation time is obtained by linear extrapolation based on the safety margin and the rate of change of the safety margin using a two-step approximation method. The warning signal is then output and the construction parameters are adjusted.

[0013] Specifically, the minimum horizontal principal stress, maximum horizontal principal stress, and vertical stress are obtained through a combination of formation testing and well logging interpretation; the original formation pore pressure is obtained through formation pressure testing; and the formation temperature is obtained through thermometry logging.

[0014] Specifically, the calibration of the tensile strength temperature decay coefficient is obtained by conducting Brazilian splitting tests on the same batch of core samples at 25℃, 150℃, 300℃, 450℃, and 600℃, respectively, and using an exponential decay model to perform least squares fitting on the temperature and rock tensile strength data, with no less than 3 samples for each temperature point.

[0015] Specifically, the calibration of the reaction kinetic parameters of the discrete activation energy distribution model is performed by thermogravimetric-mass spectrometry (TGA) experiments under multiple heating rates. The parameters are obtained by fitting the discrete activation energy distribution model with a Bayesian inference framework, which uses Markov chain Monte Carlo sampling to globally estimate the high-dimensional parameter space.

[0016] Specifically, the calibration of the strength weakening amplitude coefficient and the strength weakening rate coefficient involves immersing the core sample in a supercritical carbon dioxide environment for 1 day, 3 days, 7 days, 14 days, and 30 days. After immersion, the uniaxial compressive strength and tensile strength of the rock are measured, and the rock strength-time data are obtained by fitting the exponential decay model.

[0017] In the multi-field coupled crack initiation pressure formula, the minimum and maximum horizontal principal stresses after temperature correction are calculated using the elastic constitutive relationship through the coefficient of thermal expansion, elastic modulus, and temperature change; the tensile strength of the rock after temperature correction is calculated using the exponential decay function through the temperature decay coefficient of tensile strength; and the pyrolysis hydrocarbon generation pressurization is calculated using the pressurization conversion coefficient, rock density, hydrocarbon yield, total organic carbon content, and effective porosity.

[0018] The strength reduction coefficient is calculated by combining the strength weakening amplitude coefficient, the strength weakening rate coefficient, the heating time, the supercritical carbon dioxide injection pressure, and the volume fraction of carbonate minerals. It reflects the degree to which supercritical carbon dioxide dissolves carbonate minerals and adsorbs them onto the surface of organic matter and clay, thereby weakening the cementation strength of the rock.

[0019] The physical constraint graph attention recursive network abstracts the shale pore and fracture network into a dynamic graph structure. The graph attention layer adopts a multi-head attention mechanism, introducing physical residual jump connections within each gated loop unit time step. The numerical residual of the pyrolysis kinetic equation and the conserved residual of the continuous equation are concatenated into an additional vector and superimposed on the candidate hidden state, driving the internal iterative correction of the network until the residual is lower than the physical residual convergence threshold.

[0020] The training dataset of the physical constraint graph attention recursive network covers a sample set of parameter combinations under the conditions of total organic carbon content of 1% to 12%, carbonate mineral volume fraction of 5% to 40%, heating rate of 1℃ / min to 20℃ / min, and supercritical carbon dioxide injection pressure of 7MPa to 30MPa. It is divided into training set, validation set and test set in a ratio of 8:1:1.

[0021] The Ritz adaptive approximation algorithm is based on the Ritz variational principle of elasticity. It selects an adaptive basis function family that includes exponential thermal source term functions and Weibull distribution mineral weakening functions. Through the Rayleigh-Ritz process, it transforms the infinite-dimensional functional optimization into a finite-dimensional linear equation system. During the iteration process, it adaptively densifies or merges the basis function arrangement according to the residual distribution of the current solution.

[0022] The dynamic adjustment index function is a weighted sum of three data points: the normalized value of the safety margin, the normalized value of the rate of change of the safety margin, and the normalized value of the uncertainty of the crack initiation pressure. The weight coefficients satisfy that the sum of the three is 1. The function is determined through regression analysis of no less than 20 sets of historical construction data, with initial values ​​of 0.5, 0.3, and 0.2.

[0023] Specifically, the iterative strategy is adjusted according to the range to which the dynamic adjustment index belongs. When the dynamic adjustment index is not lower than 2, the time step is increased and an operator splitting strategy is adopted; when the dynamic adjustment index is between 1 and 2, the current setting is maintained and Newton-Krylov iteration is adopted; when the dynamic adjustment index is between 0.5 and 1, the time step is reduced and line search correction is enabled; when the dynamic adjustment index is lower than 0.5, the time step is further reduced and an early warning process is triggered.

[0024] In the two-step approximation method, the safety margin change rate is calculated by the difference between the fracture initiation pressure change rate and the bottom hole pressure change rate, and the expected fracture initiation time is obtained by linear extrapolation of the ratio of the safety margin to the absolute value of the safety margin change rate.

[0025] The warning threshold is set at 5 MPa. It is determined by subtracting twice the standard deviation from the mean of the minimum measured safety margin before the cracking event occurs, based on statistical analysis of no less than 15 sets of indoor cracking simulation test data. When applied in the field, it is adjusted after verification through no less than 3 sets of pilot tests.

[0026] The construction parameters are adjusted by reducing the supercritical carbon dioxide injection pressure or the heating rate based on the deviation between the expected crack initiation time and the warning threshold. The construction parameters refer to the collective term for four engineering parameters: supercritical carbon dioxide injection pressure, target heating temperature, heating time, and heating rate.

[0027] The system includes a readable storage medium containing program instructions that execute the aforementioned method when run on a computer.

[0028] This invention constructs a multi-field coupled fracturing pressure formula encompassing thermal field, flow field, solid mechanical field, and chemical reaction field. By combining a physical constraint graph attention recursive network, Ritz adaptive approximation algorithm, and dynamic adjustment exponential function, it achieves real-time and accurate determination of the fracturing moment in low-maturity shale under supercritical carbon dioxide heating. This solves the technical problem of being unable to determine the fracturing moment in low-maturity shale under supercritical carbon dioxide heating in real-time and accurately under multi-field coupling conditions.

[0029] The multi-field coupled fracturing pressure formula of this invention integrates the pressurization from kerogen pyrolysis, the intensity reduction from supercritical carbon dioxide, the thermal expansion stress correction, and the original geostress into a unified analytical framework, enabling quantitative characterization of the contributions of each physicochemical field. The physical constraint graph attention recursive network encodes the spatial heterogeneity of pores and fractures as graph structure features, and embeds pyrolysis kinetic constraints and fluid conservation constraints into the recursive update process through physical residual jump connections, compensating for the deficiency of the pure data-driven model in generalization under sparse sample conditions. The Ritz adaptive approximation algorithm achieves synchronous output of the spatial distribution and confidence interval of fracturing pressure through an adaptive basis function encryption strategy. The dynamically adjusted exponential function adaptively adjusts the iteration strategy and time step based on a comprehensive evaluation of the safety margin, the rate of change of the safety margin, and the uncertainty of the fracturing pressure, ensuring the computational efficiency and stability of real-time prediction.

[0030] In summary, this invention solves the technical problem mentioned in the background art of being unable to accurately determine the initiation time of fracturing in low-maturity shale under supercritical carbon dioxide heating in real time under multi-field coupling conditions. Attached Figure Description

[0031] Figure 1 This is a flowchart of the method of the present invention.

[0032] Figure 2 Curves showing the evolution of fracturing pressure and safety margin as a function of heating time in low-maturity shale under supercritical carbon dioxide heating.

[0033] Figure 3 The evolution of the confidence interval of fracturing pressure in low-maturity shale under supercritical carbon dioxide heating as a function of heating time. Detailed Implementation

[0034] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below.

[0035] like Figure 1The diagram shown is a flowchart of a method for determining the initiation of cracks in low-maturity shale heated by supercritical carbon dioxide, provided by this invention. This method includes the following steps:

[0036] S01. By integrating conventional well logging, core experiments and formation testing data, obtain the minimum horizontal principal stress, maximum horizontal principal stress, vertical stress, original formation pore pressure, formation temperature, total organic carbon content, carbonate mineral volume fraction, initial tensile strength, elastic modulus, Poisson's ratio, coefficient of thermal expansion, rock density, as well as supercritical carbon dioxide injection pressure, target heating temperature, heating time and heating rate of medium and low maturity shale formations.

[0037] S02. Based on the medium- and low-maturity shale core samples described in step S01, conduct high-temperature rock mechanics experiments, high-temperature and high-pressure closed hydrocarbon generation pressurization experiments, pyrolysis kinetics experiments, and supercritical carbon dioxide immersion mechanical comparison experiments. Calibrate the tensile strength temperature decay coefficient, pressurization conversion coefficient, discrete activation energy distribution model reaction kinetic parameters, strength weakening amplitude coefficient, and strength weakening rate coefficient, and determine the effective porosity and pore throat connectivity coefficient.

[0038] S03. Substitute the parameters described in step S01 and the calibration parameters described in step S02 into the multi-field coupled initiation pressure formula to calculate the initiation pressure at each time moment, then input the physical constraint graph attention recursive network, and output the corrected initiation pressure and safety margin.

[0039] S04. Substitute the safety margin, the rate of change of safety margin, and the uncertainty of the initiation pressure described in step S03 into the Ritz adaptive approximation algorithm to output the spatial distribution of the initiation pressure and the confidence interval of the initiation pressure at each time point.

[0040] S05. Substitute the normalized value of safety margin, the normalized value of safety margin change rate, and the normalized value of fracturing pressure uncertainty mentioned in step S04 into the dynamic adjustment index function to calculate the dynamic adjustment index. Adjust the time step and iteration strategy of the physical constraint graph attention recursive network according to the interval to which the dynamic adjustment index belongs, and update the safety margin according to the difference between the bottom hole pressure and the fracturing pressure obtained by real-time monitoring.

[0041] S06. If the safety margin described in step S05 is lower than the warning threshold, the expected crack initiation time is obtained by linear extrapolation based on the safety margin and the rate of change of the safety margin using a two-step approximation method, and a warning signal is output and the construction parameters are adjusted.

[0042] In step S01, the minimum horizontal principal stress, maximum horizontal principal stress, and vertical stress are obtained through a combination of formation testing and well logging interpretation; the original formation pore pressure is obtained through formation pressure testing; the formation temperature is obtained through thermometric logging; the total organic carbon content is obtained through combustion oxidation determination of the core sample using a carbon-sulfur analyzer; the volume fraction of carbonate minerals is obtained through X-ray diffraction mineral analysis; the initial tensile strength, elastic modulus, and Poisson's ratio of the rock are obtained through core mechanics experiments under room temperature conditions; the coefficient of thermal expansion is obtained through thermal expansion experiments of the core sample using a thermal expansion meter; the rock density is calculated by measuring the mass and volume of the core sample; the supercritical carbon dioxide injection pressure, target heating temperature, heating time, and heating rate are determined according to the engineering design scheme; the construction parameters refer to the collective term for the four engineering parameters: supercritical carbon dioxide injection pressure, target heating temperature, heating time, and heating rate.

[0043] The tensile strength temperature decay coefficient mentioned in step S02 is obtained by the following method: Brazilian splitting tests are conducted on the same batch of core samples at 25℃, 150℃, 300℃, 450℃, and 600℃ respectively, and the rock tensile strength values ​​at each temperature are recorded. The temperature and rock tensile strength data are fitted with a least squares model using an exponential decay model, with no less than 3 samples for each temperature group. The pressure conversion coefficient is obtained by the following method: Core samples with a known total organic carbon content are heated in a closed autoclave under different temperatures and time conditions. The ratio of pressure rise in the reactor to pyrolysis hydrocarbon yield was measured, normalized by rock density and effective porosity, and determined using multiple linear regression. The reaction kinetic parameters of the discrete activation energy distribution model were determined by thermogravimetric-mass spectrometry (TGA) experiments under multiple heating rates, and obtained by fitting the discrete activation energy distribution model with a Bayesian inference framework. The intensity weakening amplitude coefficient and intensity weakening rate coefficient were obtained by immersing the core sample in a supercritical carbon dioxide environment for 1, 3, 7, 14, and 30 days. After extraction, the uniaxial compressive strength and tensile strength of the rock were measured, and the rock strength-time data were obtained by fitting the exponential decay model. The effective porosity was obtained through a helium porosity measurement experiment. The pore throat connectivity coefficient was calculated by combining the pore throat diameter distribution obtained from the mercury intrusion porosimetry experiment with the specific surface area data obtained from the nitrogen adsorption experiment, reflecting the unobstructed flow of fluid between pore throats. The discrete activation energy distribution model characterizes the kerogen pyrolysis process as a superposition of a set of parallel first-order reactions with different activation energies and corresponding frequency factors, and is obtained through discretization. The activation energy distribution function describes the heterogeneity of organic matter pyrolysis; the Bayesian inference framework encodes prior geological knowledge as a prior parameter distribution, performs global estimation of the high-dimensional parameter space through Markov chain Monte Carlo sampling, and propagates the posterior parameter distribution to the fracture initiation pressure prediction to quantify the confidence interval; the Markov chain Monte Carlo sampling is a method of random walk sampling on the posterior parameter distribution by constructing a Markov chain that satisfies detailed equilibrium conditions, which is used to achieve unbiased estimation of the posterior distribution under the conditions of high parameter space dimensionality and limited experimental data.

[0044] The multi-field coupled fracture initiation pressure formula in step S03 is expressed as follows:

[0045] ;

[0046] In the formula, The initiation pressure is MPa. The reference pressure, in MPa, is the initial confining pressure set in the core mechanics experiment described in step S02. The minimum horizontal principal stress after temperature correction, in MPa The maximum horizontal principal stress after temperature correction, in MPa The tensile strength of the rock after temperature correction, in MPa This is the strength reduction factor, dimensionless. The original formation pore pressure mentioned in step S01, in MPa For the pyrolysis hydrocarbon generation pressurization, MPa; the calculation formulas for the minimum and maximum horizontal principal stresses after temperature correction are expressed as follows:

[0047] ;

[0048] In the formula, The minimum or maximum horizontal principal stress mentioned in step S01, in MPa The coefficient of thermal expansion mentioned in step S01, The elastic modulus mentioned in step S01, in MPa The temperature change, K, is calculated from the difference between the target heating temperature described in step S01 and the formation temperature described in step S01; the formula for calculating the tensile strength of the rock after temperature correction is as follows:

[0049] ;

[0050] In the formula, The initial tensile strength of the rock described in step S01, in MPa The tensile strength temperature decay coefficient mentioned in step S02, The calculation formula for the boosting of pyrolysis hydrocarbon generation is expressed as follows:

[0051] ;

[0052] In the formula, The boost conversion coefficient mentioned in step S02, The density of the rock mentioned in step S01, The hydrocarbon yield is dimensionless and is calculated from the reaction kinetic parameters of the discrete activation energy distribution model described in step S02 and the target heating temperature and heating time described in step S01. The total organic carbon content mentioned in step S01 is dimensionless. The effective porosity mentioned in step S02 is dimensionless; the formula for calculating the strength reduction factor is as follows:

[0053] ;

[0054] In the formula, The intensity weakening amplitude coefficient mentioned in step S02, The strength weakening rate coefficient mentioned in step S02, The heating time described in step S01, d The supercritical carbon dioxide injection pressure mentioned in step S01, in MPa The volume fraction of carbonate minerals mentioned in step S01 is dimensionless; the pyrolysis hydrocarbon generation pressurization refers to the increase in pore pressure caused by the accumulation of hydrocarbon gases after kerogen is pyrolyzed in a confined pore space; the strength reduction coefficient is a dimensionless coefficient reflecting the degree to which supercritical carbon dioxide weakens the cementation strength and cohesion of rocks by dissolving carbonate minerals and adsorbing onto the surface of organic matter and clay.

[0055] The specific structure of the physical constraint graph attention recursive network described in step S03 is as follows: the shale pore fracture network is abstracted into a dynamic graph structure. The feature vector of each pore throat node consists of four types of state variables: local temperature, pore pressure, carbonate mineral volume fraction, and effective porosity. The edge features consist of the pore throat radius and the pore throat connectivity coefficient. The local temperature is obtained by summing the products of the formation temperature described in step S01 and the heating rate and heating time described in step S01. The pore pressure is obtained from the intermediate results in the calculation process of the multi-field coupled fracture initiation pressure formula described in step S03. The pore throat radius is obtained through mercury intrusion porosimetry experiments. The graph attention layer adopts a multi-head attention mechanism, with each attention head calculating independently. Weighted aggregation within the node's neighborhood is used. The attention weights are obtained by normalizing the node feature vectors through a bilinear transformation and a normalized exponential function. The multi-head outputs are concatenated and then subjected to linear projection for dimensionality reduction, achieving anisotropic modeling of spatially heterogeneous mineral distribution and pore connectivity. At the graph level, node features are aggregated in parallel using global mean pooling and global max pooling to obtain the global state representation vector at the current time. This global state representation vector is input into a gated loop unit module. The gated loop unit progressively updates the hidden state along the construction time axis. The update gate and reset gate control the degree of retention of historical information and the proportion of new input information integration, respectively. Within each time step of the gated loop unit, a physical residual jump connection is introduced to connect the previous... The numerical residuals of the pyrolysis kinetic equations at each time step are concatenated with the conservation residuals of the continuity equations to form an additional vector. This additional vector is multiplied by a learnable gating weight matrix and then superimposed onto the candidate hidden states of the gated recurrent unit. This guides the hidden state update direction with physical constraints, driving iterative correction within the network until the residuals fall below the physical residual convergence threshold before exiting the inner loop. The output layer consists of two parallel fully connected sublayers that predict the initiation pressure and safety margin at the current time step. The loss function is composed of a weighted sum of the mean square error data fitting term and the fluid mass conservation residual term, with the weight coefficients determined through cross-validation. The steps for establishing the training dataset for the physical constraint graph attention recurrent network specifically include: based on the steps described in step S02... Based on experimental data and finite element numerical simulation results, a parameter combination sample set was constructed covering conditions such as total organic carbon content of 1%–12%, carbonate mineral volume fraction of 5%–40%, heating rate of 1℃ / min–20℃ / min, and supercritical carbon dioxide injection pressure of 7MPa–30MPa. The finite element numerical simulation method was used to generate time-series data of pore pressure fields and fracture initiation pressure labels under each parameter combination. Simultaneously, the experimental measurement data described in step S02 was incorporated as real samples, and the sets were divided into training, validation, and test sets in an 8:1:1 ratio. The specific steps for training the physical constraint graph attention recursive network include: using the Adam optimization algorithm, with an initial learning rate set to... The loss decays to 0.5 times its original value every 50 rounds, with a batch size of 32 and a maximum training round count of 500. An early stopping strategy monitors the validation set loss; training terminates if there is no decrease after 20 consecutive rounds. During training, the input feature vector is standardized with zero mean and unit variance, and the output layer prediction is destandardized to restore the physical dimensions. The physical residual convergence threshold is set to [value missing]. The prediction accuracy and iteration count at different threshold levels were determined through a trade-off analysis on no fewer than 10 sets of validation samples. The physical constraint graph attention recursive network encodes the spatial heterogeneity of pores and fractures as graph structure features and embeds pyrolysis kinetic constraints and fluid conservation constraints into the recursive update process using physical residual jump connections. This allows the model to make predictions under physical equation constraints even when dense monitoring data is lacking, effectively overcoming the problem of insufficient generalization ability of pure data-driven models under sparse sample conditions, and supporting real-time assessment of crack initiation pressure and safety margin at construction sites. The gated loop unit is a type of unit that updates gates and... The recurrent neural network unit structure for resetting gate control information flow is used to capture long-range dependencies in time series modeling; the global mean pooling refers to the operation of calculating the arithmetic mean of the feature vectors of all nodes in the graph to obtain the global representation; the global max pooling refers to the operation of taking the maximum value of the feature vectors of all nodes in the graph according to their dimensions to obtain the global representation; the physical residual jump connection refers to the structural design of directly introducing the physical equation residual vector into the hidden state update path of the recurrent network in a jump connection manner; the safety margin refers to the difference between the current fracturing pressure and the bottom hole pressure, reflecting the pressure margin between the current construction state and the critical fracturing state.

[0056] In step S04, the Ritz adaptive approximation algorithm is based on the Ritz variational principle of elasticity. It transforms the problem of fracture initiation pressure field distribution under multi-field coupling conditions into an energy functional extremum problem. An adaptive basis function family matching the microstructure of shale pores is selected. This family includes exponential pyrolysis source functions and Weibull-distributed mineral weakening functions. Through a Rayleigh-Ritz process, the infinite-dimensional functional optimization is transformed into a finite-dimensional linear equation system. During iteration, the basis function arrangement is refined in high-residue regions based on the residual distribution of the current solution, while sparse basis functions are merged in low-residue regions. The final output includes the spatial distribution of initiation pressure and the confidence interval of initiation pressure at each time point; the exponential pyrolysis source term function is an exponential family function constructed with the spatiotemporal distribution of the kerogen pyrolysis reaction rate as the independent variable, wherein the kerogen pyrolysis reaction rate is calculated from the reaction kinetic parameters of the discrete activation energy distribution model described in step S02 and the current temperature; the Weibull distribution-type mineral weakening function is a Weibull cumulative distribution function constructed with the spatial distribution of carbonate mineral volume fraction described in step S01 as the shape parameter, used to describe the asymmetric spatial distribution characteristics of the mineral weakening effect; the Ritz adaptive... The approximation algorithm directly encodes the physical mechanism into the construction form of an adaptive basis function family, enabling the approximation process to match the spatial distribution characteristics of the pyrolysis source term and mineral weakening. The adaptive encryption strategy concentrates computational resources on regions with large initiation pressure gradients and directly provides an estimate of the initiation pressure uncertainty through the statistical distribution of the basis function coefficients. The initiation pressure confidence interval refers to the range of values ​​within which the true initiation pressure value falls with a high probability at a given confidence level, calculated from the statistical distribution of the basis function coefficients output by the Ritz adaptive approximation algorithm. The Ritz variational principle of elasticity refers to... The variational method for solving elasticity problems involves taking the extreme values ​​of the potential energy functional within the displacement field function space of the boundary conditions. The Rayleigh-Ritz process refers to representing the approximate solution as a linear combination of a set of basis functions, and transforming the infinite-dimensional variational problem into a finite-dimensional linear algebraic equation system by setting the partial derivatives of the functional with respect to the coefficients of each basis function to zero. The uncertainty of the initiation pressure refers to the statistical standard deviation of the pressure values ​​at each point in the spatial distribution of the initiation pressure output by the Ritz adaptive approximation algorithm, reflecting the confidence level of the model prediction results. After being output in step S04, it is passed to step S05 for calculating the dynamic adjustment index.

[0057] The calculation formula for the dynamic adjustment exponential function described in step S05 is as follows:

[0058] ;

[0059] In the formula, It is a dynamic adjustment index, dimensionless. For safety margin, MPa For reference safety margin, MPa is taken as the safety margin calculated at the start of construction in step S03. For the rate of change of safety margin, For reference, the rate of change of safety margin The value is the root mean square of the rate of change of the safety margin over the first 10 time steps. The uncertainty of the initiation pressure mentioned in step S04, in MPa For reference, the initiation pressure uncertainty, MPa, is taken as the initiation pressure uncertainty output at the initial time in step S04. , , The weighting coefficients are dimensionless and satisfy the following conditions: The initial value was determined through regression analysis of no fewer than 20 sets of historical construction data. , , ;when At that time, the time step is increased by a factor of 1.5, and the iteration convergence threshold is set to... The operator splitting strategy is adopted; when When the current time step and iteration convergence threshold are maintained, Newton-Krylov iteration is used; when When the time step is reduced to 0.5 times the original time step, the iteration convergence threshold is set to... Enable line search correction; when When the time step is reduced to 0.25 times the original time step, the iteration convergence threshold is set to... Simultaneously triggering step S06, the early warning process; the dynamic adjustment exponential function is a weighted summation function of three data items: the normalized value of the safety margin, the normalized value of the rate of change of the safety margin, and the normalized value of the uncertainty of the crack initiation pressure, used to quantify the degree to which the current construction state deviates from the safe working condition; the operator splitting strategy refers to decomposing the four-field coupling equation system of heat-fluid-solid-chemical into several sub-problems according to the physical field, and solving each sub-problem sequentially in each time step, with the output of the previous sub-problem serving as the known input of the next sub-problem, and approximating the fully coupled solution through multiple rounds of sub-step iteration; the Newton-Krylov iteration replaces the solution of the linear equation system of Newton iteration with an iterative linear solver based on the Krylov subspace, avoiding explicit assembly and direct solution of large-scale Jacobian matrices, and is suitable for multi-field coupling equation systems caused by cross terms. The solution method for asymmetric sparse Jacobian structures; the line search correction refers to performing a step size search in each update direction of the Newton-Krylov iteration to ensure that the objective function value decreases after each iteration, which is a correction strategy to prevent iteration divergence when the safety margin change rate is large; the first 10 time steps refer to the first 10 calculation time steps from the start of construction, and the reference safety margin change rate is calculated by taking the root mean square of the safety margin change rate sequence of the first 10 time steps; the safety margin normalized value refers to the ratio of the safety margin to the reference safety margin; the safety margin change rate normalized value refers to the ratio of the absolute value of the safety margin change rate to the reference safety margin change rate; the crack initiation pressure uncertainty normalized value refers to the ratio of the crack initiation pressure uncertainty in step S04 to the reference crack initiation pressure uncertainty.

[0060] In the two-step approximation method described in step S06, the formula for calculating the rate of change of safety margin is as follows:

[0061] ;

[0062] The formula for calculating the expected crack initiation time is expressed as follows:

[0063] ;

[0064] In the formula, To predict the crack initiation time, d At the current moment, d For reference time scale, d is the heating time described in step S01. The bottom hole pressure obtained through real-time monitoring as described in step S05, in MPa This represents the rate of change of crack initiation pressure. The fracture initiation pressure difference between adjacent time moments is calculated using the multi-field coupled fracture initiation pressure formula described in step S03. This represents the rate of change of bottom hole pressure. The warning threshold is calculated from the pressure difference at adjacent wellbore times obtained by real-time monitoring; the warning threshold is 5 MPa, determined by subtracting twice the standard deviation from the mean of the minimum measured safety margin before the occurrence of the fracturing event through statistical analysis of no less than 15 sets of indoor fracturing simulation experimental data, and adjusted after verification by no less than 3 sets of pilot tests in field application; the construction parameters are adjusted by reducing the supercritical carbon dioxide injection pressure or reducing the heating rate based on the deviation between the expected fracturing time and the warning threshold.

[0065] Optionally, the present invention also provides a computer-based method for forming a supercritical carbon dioxide heating crack initiation judgment system for low-maturity shale, wherein the computer is equipped with a readable storage medium storing program instructions, and the program instructions can execute the above-described method when running on the computer.

[0066] The specific implementation of step S01 is as follows: First, using conventional logging curves such as sonic logging, density logging, and resistivity logging, combined with formation test data, the minimum horizontal principal stress, maximum horizontal principal stress, and vertical stress of the formation are obtained through joint inversion using logging interpretation methods. The vertical stress is calculated by integrating the density of the overlying formation, and the horizontal principal stress is constrained by a combination of fracturing tests and logging elastic parameters. The original formation pore pressure is directly measured in the target section using a formation pressure testing instrument while drilling. The formation temperature is determined using cable temperature logging under stable heat flow conditions in the wellbore. The total organic carbon content is obtained by measuring the carbon content after complete combustion and oxidation of the core sample using a carbon-sulfur analyzer. The volume fraction of carbonate minerals is quantitatively determined by X-ray diffraction mineral analysis. The initial tensile strength, elastic modulus, and Poisson's ratio of the rock are determined by uniaxial or triaxial mechanical experiments at room temperature. The coefficient of thermal expansion is calculated by continuously measuring the change in the length of the core sample during heating using a thermal dilatometer. The rock density is calculated by the ratio of the dried mass of the core sample to its geometric volume. The four construction parameters—supercritical carbon dioxide injection pressure, target heating temperature, heating time, and heating rate—were determined based on the engineering design scheme. These parameters constitute the entire input basis for subsequent multi-field coupling calculations.

[0067] The specific implementation of step S02 involves conducting four types of calibration experiments based on the core samples obtained in step S01. The high-temperature rock mechanics experiment uses the Brazilian splitting method, taking at least three samples from each of five temperature points (25℃, 150℃, 300℃, 450℃, and 600℃) from the same batch of core samples to determine tensile strength. The temperature-strength data are then fitted using an exponential decay model to obtain the temperature decay coefficient of tensile strength. The high-temperature, high-pressure closed-loop hydrocarbon generation pressurization experiment involves heating core samples with a known total organic carbon content in a sealed autoclave under multiple temperature and time conditions. The pressure rise and pyrolysis hydrocarbon yield are recorded. After normalization using rock density and effective porosity, the pressurization conversion coefficient is determined through multiple linear regression. Thermogravimetric-mass spectrometry (TGA-MS) was used to analyze the organic matter in the core samples under multiple heating rates. The pyrolysis process of kerogen was characterized as a superposition of parallel first-order reactions with different activation energies. Using a Bayesian inference framework and Markov chain Monte Carlo sampling, the high-dimensional parameter space was globally estimated, and the reaction kinetic parameters and their posterior distributions of the discrete activation energy distribution model were fitted. A mechanical comparison experiment after supercritical carbon dioxide immersion was conducted by placing the core samples in a supercritical environment. In the environment, uniaxial compressive and tensile strengths were measured at five time points: 1 day, 3 days, 7 days, 14 days, and 30 days. The strength-time curves were fitted using an exponential decay model to obtain the strength weakening amplitude coefficient and the strength weakening rate coefficient. Effective porosity was obtained through helium porosimetry. The pore-throat connectivity coefficient was calculated by combining the pore-throat diameter distribution obtained from mercury intrusion porosimetry with the specific surface area data obtained from nitrogen adsorption experiments, reflecting the unobstructed flow of fluid between pores and throats. These calibration parameters provide quantitative physical constraints for the multi-field coupling calculations in step S03.

[0068] The specific implementation of step S03 involves substituting all parameters obtained in steps S01 and S02 into the multi-field coupled fracturing pressure formula and calculating the fracturing pressure time-by-time. The multi-field coupled fracturing pressure formula is based on the tensile failure criterion of elasticity. It corrects the minimum and maximum horizontal principal stresses through the product of the elastic constitutive relation, thermal expansion coefficient, elastic modulus, and temperature change. It calculates the temperature-corrected rock tensile strength using an exponential decay function and a tensile strength temperature decay coefficient. It calculates the hydrocarbon yield using reaction kinetic parameters from a discrete activation energy distribution model, the target heating temperature, and heating time, and then calculates the pyrolysis hydrocarbon generation pressurization by combining the pressurization conversion coefficient, rock density, total organic carbon content, and effective porosity. Finally, it calculates the strength reduction coefficient using the strength weakening amplitude coefficient, strength weakening rate coefficient, heating time, supercritical carbon dioxide injection pressure, and carbonate mineral volume fraction. The calculated fracturing pressures at each time point are input into a physical constraint graph attention recursive network. This network abstracts the shale pore-fracture network into a dynamic graph structure. Each pore throat node's feature vector consists of local temperature, pore pressure, carbonate mineral volume fraction, and effective porosity, while edge features consist of pore throat radius and pore throat connectivity coefficient. The graph attention layer employs a multi-head attention mechanism to anisotropically model the spatially heterogeneous mineral distribution and pore connectivity. Global mean pooling and global max pooling are aggregated in parallel and then input into a gated recurrent unit module to progressively update the hidden state along the time axis. At each time step, a physical residual jump connection is introduced, concatenating the numerical residual of the pyrolysis kinetic equation with the conservation residual of the continuity equation into an additional vector, which is then superimposed onto the candidate hidden state, driving iterative correction within the network until the residual falls below the physical residual convergence threshold. The output layer outputs the corrected fracture initiation pressure and safety margin at the current moment in parallel. The Adam optimization algorithm is used during training, with an initial learning rate of... The loss decreases to 0.5 times every 50 rounds, the batch size is 32, the maximum number of training rounds is 500, the early stopping strategy monitors the loss of the validation set, and the program terminates if there is no decrease for 20 consecutive rounds.

[0069] The specific implementation of step S04 involves substituting the safety margin, rate of change of safety margin, and uncertainty of initiation pressure output from step S03 into the Ritz adaptive approximation algorithm. Based on the Ritz variational principle of elasticity, the problem of initiation pressure field distribution under multi-field coupling conditions is transformed into an energy functional extremum problem. An adaptive basis function family matching the microstructure of shale pores is selected. This family includes an exponential pyrolysis source term function with the spatiotemporal distribution of kerogen pyrolysis reaction rate as the independent variable, and a Weibull distribution-type mineral weakening function with the spatial distribution of carbonate mineral volume fraction as the shape parameter. The infinite-dimensional functional optimization is transformed into a finite-dimensional linear equation system through a Rayleigh-Ritz process. During iteration, the basis function arrangement is refined in high-residue regions based on the residual distribution of the current solution, while sparse basis functions are merged in low-residue regions, concentrating computational resources on regions with large initiation pressure gradients. Finally, the spatial distribution of initiation pressure and the confidence interval of initiation pressure at each time point are output, where the uncertainty of initiation pressure is directly provided by the statistical standard deviation of the basis function coefficients.

[0070] The specific implementation of step S05 is as follows: the normalized value of the safety margin, the normalized value of the rate of change of the safety margin, and the normalized value of the uncertainty of the initiation pressure output in step S04 are substituted into the dynamic adjustment index function. The three normalized values ​​are normalized with reference quantities such as the safety margin at the start of construction, the root mean square of the rate of change of the safety margin in the first 10 time steps, and the uncertainty of the initiation pressure at the initial time. The initial values ​​of the weighting coefficients are 0.5, 0.3, and 0.2, and are updated through regression of historical data. According to the range of the dynamic adjustment index, when the dynamic adjustment index is not less than 2, the time step is increased by 1.5 times and an operator splitting strategy is adopted to decompose the four-field equation system of heat-fluid-solid-chemical into sequential subproblems and solve them sequentially; when the dynamic adjustment index is between 1 and 2, the current setting is maintained and Newton-Krylov iteration is adopted; when the dynamic adjustment index is between 0.5 and 1, the time step is reduced to 0.5 times and line search correction is enabled; when the dynamic adjustment index is less than 0.5, the time step is further reduced to 0.25 times and the early warning process of step S06 is triggered. At the same time, the safety margin is updated based on the difference between the bottom hole pressure and the fracturing pressure obtained from real-time monitoring.

[0071] The specific implementation of step S06 is as follows: when the safety margin is lower than the warning threshold of 5 MPa, a two-step approximation method is used to linearly extrapolate the expected fracturing time. The first step calculates the safety margin change rate by obtaining the rate of change of the safety margin at the current moment through the difference between the rate of change of the fracturing pressure and the rate of change of the bottom hole pressure monitored in real time. The second step divides the current safety margin by the absolute value of the safety margin change rate, normalizes it using the heating time as a reference time scale, and obtains the expected fracturing time. After the expected fracturing time is determined, the system outputs a warning signal. Based on the deviation between the expected fracturing time and the warning threshold, the construction personnel adjust the construction parameters by reducing the supercritical carbon dioxide injection pressure or the heating rate. The warning threshold of 5 MPa is determined by statistical analysis of no less than 15 sets of indoor fracturing simulation experimental data, by subtracting twice the standard deviation from the mean of the minimum measured safety margin before the fracturing event. In field application, it is adjusted after verification through no less than 3 sets of pilot tests.

[0072] It should be noted that the key technologies of this invention include: the multi-field coupled fracturing pressure formula incorporates pyrolysis hydrocarbon generation pressurization, supercritical carbon dioxide intensity reduction, thermal expansion stress correction, and original geostress into a unified analytical framework, enabling quantitative coupling of the contributions of each physicochemical field in the same calculation system, significantly improving the ability to describe complex multi-field evolution processes compared to single-field or dual-field models; the physical constraint graph attention recursive network embeds the physical equation residuals into the hidden state update path of the gated recurrent unit in a skip connection manner, enabling the model to make effective predictions based on physical constraints even under sparse sample conditions, overcoming the inherent defect of the pure data-driven model having severely degraded generalization ability when training data is insufficient; and the Ritz adaptive approximation algorithm, through adaptive basis function encryption and sparse merging strategies, simultaneously outputs the spatial distribution of fracturing pressure and uncertainty estimates, providing a reliable statistical basis for the calculation of the dynamic adjustment index. Three key technologies work together: the multi-field coupling formula provides physically consistent input and residual constraints for the physical constraint graph attention recursive network, ensuring that the corrected crack initiation pressure output by the network satisfies multi-field physical self-consistency; the safety margin output by the physical constraint graph attention recursive network and the uncertainty output by the Ritz adaptive approximation algorithm together constitute the computational input for the dynamic adjustment index, enabling the adaptive control of the iterative strategy to simultaneously consider prediction accuracy and computational stability; the synergistic effect of these three technologies ensures that the entire judgment process maintains reliable crack initiation prediction capability under the triple constraints of strong multi-field coupling, sparse samples, and real-time computation.

[0073] It should be noted that this invention also solves the following technical problem: In the construction process of supercritical carbon dioxide heating of low-maturity shale, due to the inherent experimental measurement uncertainties of key parameters such as pyrolysis kinetic parameters and mineral weakening coefficients, existing methods cannot provide quantitative confidence interval estimates for the initiation pressure prediction results. This makes it difficult for construction personnel to assess the reliability of the prediction results, and consequently, to formulate reasonable safety margins for construction. This invention introduces a Bayesian inference framework into the parameter calibration process of a discrete activation energy distribution model, obtains the posterior distribution of key parameters using Markov chain Monte Carlo sampling, and propagates the parameter uncertainties to the initiation pressure prediction through a multi-field coupling formula. Then, the statistical distribution of the basis function coefficients of the Ritz adaptive approximation algorithm directly provides estimates of the initiation pressure uncertainty at each time step, thereby outputting a confidence interval for the initiation pressure with a clear confidence level. This solves the technical problem that existing methods cannot quantitatively characterize the uncertainty of initiation pressure prediction, and provides a quantifiable probabilistic basis for the reasonable setting of construction safety margins.

[0074] Specifically, the principle of this invention is as follows: The fundamental reason why this invention can solve the above-mentioned technical problems lies in its technical solution, which starts from the physical mechanism, quantifies the multi-field coupling effect step by step, and uses artificial intelligence and variational approximation methods to collaboratively handle model uncertainties and real-time calculation requirements. The multi-field coupling initiation pressure formula is based on the elastic mechanical fracture criterion. It converts the hydrocarbon yield and pressure conversion coefficient obtained by calculating the pyrolysis hydrocarbon generation pressure through the discrete activation energy distribution model kinetic equation into the pore pressure increment. It quantitatively expresses the supercritical carbon dioxide intensity reduction through the carbonate mineral dissolution and organic matter adsorption mechanism using an exponential decay function. It corrects the geostress through the elastic constitutive relation of thermal expansion stress, thereby ensuring that the formula is physically self-consistent. The physical constraint graph attention recursive network guides the direction of model parameter update by using the physical equation residual as a regularization constraint under data sparsity conditions, so that the prediction results always satisfy the mass conservation and pyrolysis kinetic constraints, avoiding the overfitting or underfitting problems caused by insufficient training samples in the pure data-driven model. The Ritz adaptive approximation algorithm transforms the problem of crack initiation pressure field distribution into an energy functional extremum problem. By constructing an adaptive basis function family that matches the physical mechanism, it directly provides uncertainty estimates, offering reliable input for the calculation of the dynamic adjustment exponent. The dynamic adjustment exponent function comprehensively quantifies the safety margin, the rate of change of the safety margin, and the uncertainty into a single exponent, adaptively adjusting the iterative strategy to ensure a dynamic balance between prediction accuracy and computational efficiency throughout the entire construction process, thereby achieving real-time and accurate determination of the crack initiation moment.

[0075] The following provides a specific embodiment 1 of the present invention, and the specific implementation of each step in this embodiment 1 is described in detail below.

[0076] The specific implementation of step S01 is to obtain the minimum horizontal principal stress through a combination of formation testing and well logging interpretation. Maximum horizontal principal stress With vertical stress All units are MPa; original formation pore pressure (MPa) is measured directly using a drilling formation pressure testing instrument; formation temperature (°C) obtained through cable temperature logging; Total organic carbon content (Dimensionless) Determined by combustion oxidation using a carbon-sulfur analyzer; volume fraction of carbonate minerals. (Dimensionless) Quantitative determination of initial tensile strength of rock by X-ray diffraction mineral analysis. (MPa), elastic modulus (MPa) and Poisson's ratio (Dimensionless) Determined by triaxial or uniaxial mechanical experiments at room temperature; coefficient of thermal expansion ( The rock density was calculated by continuously measuring the length change of the core sample during the heating process using a thermal dilatometer. ( ) Calculated by the ratio of mass to geometric volume of dried core sample; supercritical Injection pressure (MPa), target heating temperature (°C), heating time (d) with heating rate (℃ / min) is determined based on the engineering design scheme, representing the temperature change. (K) is calculated by the following formula:

[0077] ;

[0078] In the formula, The target heating temperature (°C) is specified. The original stratum temperature (°C) is given. The difference between the two is the temperature rise experienced by the formation during the heating process. Here, the difference between °C and K is the same. Numerically, it is equal to the difference in degrees Celsius between the two temperatures.

[0079] The specific implementation of step S02 is the tensile strength temperature decay coefficient. ( The tensile strength of the rock was obtained through the Brazilian splitting test. At least three samples were taken from each of five temperature points (25℃, 150℃, 300℃, 450℃, and 600℃) of the same batch of core samples, and the tensile strength of the rock at each temperature was measured. (MPa), the least squares fit is performed using an exponential decay model, and the fitting equation is expressed as follows:

[0080] ;

[0081] In the formula, For the first The tensile strength (MPa) of the rock measured at several temperature points. The initial tensile strength of the rock at room temperature (MPa). For the first Each experimental temperature (K) The reference temperature is room temperature, 298K. The temperature decay coefficient of the tensile strength to be fitted ( The objective function for the sum of squared residuals of all temperature data points is obtained by minimizing the objective function, which is expressed as follows:

[0082] ;

[0083] In the formula, The total number of experimental samples (dimensionless). The objective function is the sum of squared residuals (dimensionless). This is achieved by... Find the derivative and set it to zero, or use numerical optimization methods to make it... Minimize the temperature decay coefficient of the tensile strength to obtain the optimal value. Boost conversion coefficient ( This was obtained through experiments in a closed autoclave, under different temperature and time conditions, for known... Core heating was performed, and the pressure rise inside the reactor was recorded. (MPa) and pyrolysis hydrocarbon yield (Dimensionless), combined with rock density With effective porosity (Dimensionless) normalization, determined by multiple linear regression. The regression equation is expressed as follows:

[0084] ;

[0085] In the formula, For the first Pressure rise (MPa) inside the experimental vessel. The reference pressure (MPa) is the initial confining pressure set in the core mechanics experiment. For the first The yield of pyrolysis hydrocarbons corresponding to each group of experiments (dimensionless). The regression error term (dimensionless) reflects the experimental measurement bias and model assumption error. By minimizing The mean square error was determined. The reaction kinetic parameters of the discrete activation energy distribution model were determined by thermogravimetric-mass spectrometry (TGA) experiments under multiple heating rates. The hydrocarbon yield from kerogen pyrolysis was also determined. (Dimensionless) Characterized by a superposition of a set of parallel first-order reactions, the calculation formula is expressed as follows:

[0086] ;

[0087] In the formula, The number of discrete activation energies (dimensionless). For the first The mass weights (dimensionless) of the group reactions satisfy... For the first Frequency factor of group response ( ) For the first Activation energy of group reaction ( ) For an ideal gas constant, the value is 8.314. Dumb variables for integration The formation temperature (K) at that time. From the start of heating to the current moment The time dummy variable (s) in the process is used to describe the cumulative effect of the pyrolysis reaction throughout the heating history; the above parameter set The distribution was obtained by fitting a Bayesian inference framework combined with Markov chain Monte Carlo sampling. The prior distribution was encoded by geological experience, and the posterior distribution was estimated through Markov chain Monte Carlo random walk sampling. The sampling process satisfied detailed equilibrium conditions. (Intensity weakening amplitude coefficient) ( ) and intensity weakening rate coefficient ( ) through supercritical The core samples were obtained by immersion experiments, in which the core samples were immersed in supercritical fluid. In the environment, the tensile strength of the rock was measured after it was taken out at five time points: 1 day, 3 days, 7 days, 14 days, and 30 days. (MPa), the intensity-time data were fitted using an exponential decay model, and the fitting equation is expressed as follows:

[0088] ;

[0089] In the formula, For the first The tensile strength (MPa) of the rock was measured after it was removed from the immersion site at each soaking time point. For the first Each soaking time point (d). and The parameters to be fitted are obtained through nonlinear least-squares fitting of intensity data at five time points. This equation describes the supercritical... The weakening process of dissolving carbonate minerals and adsorbing them onto organic matter surfaces over time, an exponential term. The time-saturation characteristics characterizing the weakening effect. Effective porosity. Pore ​​connectivity coefficient was obtained through helium gas porosity measurement experiments; (Dimensionless) The diameter distribution of the pore throat was obtained through mercury intrusion porosimetry, and the specific surface area was obtained through nitrogen adsorption experiments. ( The calculation is performed jointly, and the formula is expressed as follows:

[0090] ;

[0091] In the formula, For the first Diameter of the larynx (m). For the corresponding frequency of occurrence (dimensionless). Normalization coefficient ( Its dimensions are in the order of molecules. Dimensions Divide by the denominator middle dimension and The product of dimensions is determined, so that Maintain dimensionless, Determined by normalization fitting of experimental data.

[0092] The specific implementation of step S03 is to substitute all the parameters obtained in steps S01 and S02 into the multi-field coupled fracturing pressure formula, which is expressed as follows:

[0093] ;

[0094] In the formula, The crack initiation pressure (MPa). This represents the minimum horizontal principal stress (MPa) after temperature correction. This represents the maximum horizontal principal stress (MPa) after temperature correction. This represents the temperature-corrected tensile strength of the rock (MPa). This is the strength reduction factor (dimensionless). The pressure increase (MPa) for pyrolysis hydrocarbon generation; this formula is derived from the tensile failure criterion of elasticity. The difference between the minimum and maximum principal stresses (three times the minimum principal stress) in the numerator on the right side of the formula reflects the stress effect at the far-field location. The tensile strength term, corrected by the strength reduction factor, reflects the supercritical stress. Weakening effect and The sum reflects the dual sources of pore pressure, with each term having the dimension of MPa. Dividing the whole by... Both sides are dimensionless. The formula for calculating the horizontal principal stress after temperature correction is as follows:

[0095] ;

[0096] In the formula, The initial value of the minimum or maximum horizontal principal stress (MPa). The second term on the right represents the horizontal principal stress (MPa) after temperature correction. Dimensions are , divided by (MPa) is followed by a dimensionless value, which is consistent with the left side. Dimensionally consistent; this formula is derived based on the constitutive relation of linear elastic thermal expansion, and the thermal expansion stress increment This reflects the superposition effect of constrained thermal stress caused by temperature rise. The formula for calculating the tensile strength of rock after temperature correction is as follows:

[0097] ;

[0098] In the formula, The temperature-corrected tensile strength of the rock (MPa) is represented by the exponent term on the right. Dimensions are , is dimensionless, and is on the left side The dimensions are consistent. The calculation formula for the pressurization of pyrolysis hydrocarbon generation is expressed as follows:

[0099] ;

[0100] In the formula, the molecule Dimensions are ( and (dimensionless), denominator The dimensionless measurement is MPa ( (dimensionless), the ratio of the two is dimensionless, and the left side The dimensions are consistent. The formula for calculating the strength reduction factor is as follows:

[0101] ;

[0102] In the formula, ( )and Multiplying the product of (dimensionless, dimensionless, MPa / MPa) by (dimensionless) and (Dimensionless) The dimensionless dimension is (dimensionless), and the whole Dimensionless; this formula describes supercritical... The weakening process of dissolving carbonate minerals and adsorbing them onto organic surfaces over time. In a physical constraint graph attention recursive network, each pore throat node... eigenvectors The statement is as follows:

[0103] ;

[0104] In the formula, For nodes The local temperature (K) is determined by The calculation yielded, where The heating rate is expressed in °C / min. The heating time (d) is here. After conversion to K and Add For nodes The pore pressure (MPa) was obtained from intermediate results of the multi-field coupling formula. For nodes Volume fraction of carbonate minerals (dimensionless) For nodes Effective porosity (dimensionless); superscript This represents the transpose of a vector. Edge eigenvectors. The statement is as follows:

[0105] ;

[0106] In the formula, For nodes With nodes The corresponding throat radius (m) between them. For nodes With nodes The dimensionless pore throat connectivity coefficients (corresponding to the pore throats) are taken from the pore throat connectivity coefficients measured in step S02. Local assignment on the corresponding aperture throat pair; multi-head attention weights The formula for calculating (dimensionless) is as follows:

[0107] ;

[0108] In the formula, For the first Bilinear weight matrix (dimensionless) for each attention head. For nodes The set of neighboring nodes, with the denominator being the pairs of nodes. The summation of the normalized exponential function of all neighboring nodes makes Satisfy normalization conditions In physical residual skip connections, the candidate hidden state update formula is expressed as follows:

[0109] ;

[0110] In the formula, for Candidate hidden state vector at time step. for The gate output vector is reset at any time (dimensionless). This is the hidden state vector from the previous time step. The input vector at the current time step, For element-wise multiplication, The candidate hidden state weight matrix is... For learnable gated weight matrix, For bias vectors, The physical residual vector from the previous time step is formed by concatenating the numerical residuals of the pyrolysis kinetic equations and the conservation residuals of the continuity equations. This process is repeated within the network until... Exit the inner loop, where Let L2 represent the norm of a vector. The loss function is expressed as follows:

[0111] ;

[0112] In the formula, Total loss (dimensionless). and The weighting coefficients (dimensionless) for the two losses satisfy the following conditions: The empirical initial value is determined on the training set through cross-validation. , The mean square error data fitting term (dimensionless) is expressed as follows:

[0113] ;

[0114] In the formula, This represents the batch sample size (dimensionless). Predict the initiation pressure (MPa) for the network. The difference between the finite element simulation or experimental label crack initiation pressure (MPa) and the pressure is divided by... The latter is dimensionless. The fluid mass conservation residual term (dimensionless) is expressed as follows:

[0115] ;

[0116] In the formula, For the first The physical residual vector of each sample (dimensionless). For reference residual scale (dimensionless), empirical values ​​are taken as follows: , It is dimensionless.

[0117] The specific implementation of step S04 is that the Ritz adaptive approximation algorithm distributes the crack initiation pressure field. (MPa) represents the family of adaptive basis functions. A linear combination of , expressed as follows:

[0118] ;

[0119] In the formula, Let m be the spatial coordinate vector. For the first Dimensional coefficients (MPa) of each basis function. The total number of basis functions (dimensionless). For dimensionless basis functions, the left side With the right side All dimensions are dimensionless; by making the energy functional pair The partial derivatives are zero, resulting in a finite-dimensional linear system of equations. The elements of the stiffness matrix and the elements of the load vector are described as follows:

[0120] ;

[0121] ;

[0122] In the formula, The integration domain is the target shale formation region ( ), The target value (MPa) of the initiation pressure at each spatial point is output by the attention recursive network of the physical constraint graph in step S03. The stiffness matrix is ​​the first Line number The column elements have the following dimensions: Dimensionless after normalization (after integration) Dimensions are ,molecular Dimensionless, denominator dimension The overall dimension is , correspond The dimensions must be consistent. Dimensions are coefficient vector The dimensionless quantity is MPa, so that Dimensions are ,and (Consistent); the adaptive basis function family includes exponential pyrolysis source term functions Compared with the weakening function of Weibull distribution type minerals The exponential pyrolysis source term function is expressed as follows:

[0123] ;

[0124] In the formula, For the first Activation energy of group reaction ( ), For the corresponding frequency factor ( ), For position Place The formation temperature (K) at that time. Reference frequency factor ( Take each group The mean of, makes It is dimensionless; the weakening function of Weibull distribution type minerals is expressed as follows:

[0125] ;

[0126] In the formula, is a scale parameter (dimensionless). For shape parameters (dimensionless). For position Volume fraction of carbonate minerals (dimensionless), overall It is dimensionless; during the iteration process, the basis functions are densified in the high residual region and merged in the low residual region according to the residual distribution of the current solution; the uncertainty of the crack initiation pressure. (MPa) is calculated from the posterior statistical standard deviation of the basis function coefficients through error propagation, and is expressed as follows:

[0127] ;

[0128] In the formula, For the first basis function coefficients The posterior standard deviation (MPa). These are dimensionless basis function values. Since it is dimensionless, all terms within the square root of the whole and their sum are also dimensionless, as shown on the left. The dimensions are consistent.

[0129] The specific implementation of step S05 is as follows: the dynamic adjustment exponential function is expressed as follows:

[0130] ;

[0131] In the formula, It is a dynamic adjustment index (dimensionless). For safety margin (MPa). For reference safety margin (MPa) For reference, the rate of change of safety margin (MPa / d) Weighting factor for reference initiation pressure uncertainty (MPa) , , satisfy The initial value is taken as , , The data was updated through regression analysis of no fewer than 20 sets of historical construction data; the three items were the ratio of the safety margin to the reference safety margin, the absolute value of the rate of change of the safety margin, and so on. The ratio of (MPa / d) to the rate of change of the reference safety margin, and the ratio of the uncertainty of the initiation pressure to the uncertainty of the reference initiation pressure, are both dimensionless; based on The time step and iteration strategy are adaptively adjusted within the specified interval, and the safety margin is determined by real-time monitoring of the bottom hole pressure. The difference between (MPa) and crack initiation pressure is continuously updated.

[0132] The specific implementation of step S06 is as follows: the normalized calculation formula for the rate of change of safety margin is expressed as follows:

[0133] ;

[0134] In the formula, The rate of change of initiation pressure (MPa / d) is calculated from the difference in initiation pressure between adjacent time points. The numerator represents the rate of change of bottom hole pressure (MPa / d), calculated from the pressure difference between adjacent time points using real-time monitoring. The dimension is MPa / d, divided by the rate of change of the reference safety margin. (MPa / d) After that, the overall value is dimensionless; the formula for calculating the expected crack initiation time is expressed as follows:

[0135] ;

[0136] In the formula, The expected crack initiation time (d) is given. For the current time (d), The heating time is taken as the reference time scale (d). , The current safety margin (MPa) is shown on the left. It is dimensionless, and the molecule on the right side is... Dimensionless, denominator It is dimensionless, with consistent dimensions on both sides; the default warning threshold is 5 MPa. When the safety margin is lower than this threshold, a warning signal is output, and the supercritical pressure is reduced based on the deviation between the expected crack initiation time and the warning threshold. Adjust the construction parameters by injecting pressure or reducing the heating rate.

[0137] To better understand and implement this invention, the following is a specific application scenario of the invention, Example 2: To illustrate the application process of the technical solution of this invention, technicians conducted a crack initiation judgment application test in a medium-low maturity continental shale stratum. The target stratum was buried at a depth of about 2800m, the stratum temperature was about 85℃, the total organic carbon content was about 6.3%, the volume fraction of carbonate minerals was about 22%, and the effective porosity was about 4.1%.

[0138] First, step S01 is executed, obtaining formation mechanical parameters through a combination of well logging interpretation and formation testing. The minimum horizontal principal stress is approximately 43 MPa, the maximum horizontal principal stress is approximately 58 MPa, the vertical stress is approximately 67 MPa, and the original formation pore pressure is approximately 28 MPa. The engineering design scheme is set to supercritical. The injection pressure was 18 MPa, the target heating temperature was 350℃, the heating time was 30 days, and the heating rate was 5℃ / min. The initial tensile strength of the rock was approximately 7.2 MPa, the elastic modulus was approximately 21 GPa, the Poisson's ratio was approximately 0.24, and the coefficient of thermal expansion was approximately... The rock has a density of approximately 2.48. The pore-throat connectivity coefficient is approximately 0.63.

[0139] Subsequently, step S02 was performed to conduct four types of calibration experiments. Three samples from each of the same batch of core samples were taken at 25℃, 100℃, 200℃, 300℃, and 400℃ to conduct Brazilian splitting tests. The temperature decay coefficient of tensile strength was obtained by fitting. The closed-loop hydrocarbon generation pressurization experiment determined the pressurization conversion coefficient to be approximately [missing value] using multiple linear regression. MPa The pyrolysis kinetics experiment yielded reaction kinetic parameters for a discrete activation energy distribution model by fitting a Bayesian inference framework. The main peak activation energy was approximately 214 kJ / m². The corresponding frequency factor is approximately Supercritical The immersion experiment was conducted at five time points, and the strength weakening amplitude coefficient obtained by fitting was approximately 0.18. The rate of weakening intensity is approximately 0.09. The calibration results for each experiment are shown in Table 1.

[0140] Table 1 Summary of key calibration parameters for step S02

[0141]

[0142] In step S03, the above parameters are substituted into the multi-field coupled fracturing pressure formula. Taking day 10 as an example: after thermal expansion stress correction, the minimum horizontal principal stress is approximately 48.7 MPa, and the maximum horizontal principal stress is approximately 65.5 MPa; after temperature correction, the tensile strength is approximately 4.8 MPa; calculated by the discrete activation energy distribution model, the hydrocarbon yield on day 10 is approximately 0.041, and the pyrolysis hydrocarbon generation pressure increase is approximately 3.2 MPa; the strength reduction factor is approximately 0.87; the comprehensive calculation yields a fracturing pressure on day 10 of approximately 39.6 MPa, with a safety margin of approximately 21.4 MPa. The evolution trend of fracturing pressure at each time point with heating time is as follows: Figure 2 As shown in the figure, the nonlinear change process of the initiation pressure continuously decreasing over time is illustrated. The deviation between the initiation pressure after correction by the physical constraint graph attention recursive network and the direct calculation result of the formula is significantly reduced in the later stages of construction, reflecting the network's effective ability to capture multi-field coupled nonlinear effects.

[0143] In step S04, the Ritz adaptive approximation algorithm, based on the Ritz variational principle of elasticity, incorporates the exponential pyrolysis source term function and the Weibull distribution mineral weakening function into an adaptive basis function family, iteratively refining the basis function arrangement in the high residual region. On day 10, the output fracture initiation pressure uncertainty is approximately 1.7 MPa, with a 95% confidence interval of 36.3 MPa to 43.0 MPa. The evolution of the fracture initiation pressure confidence interval with heating time is as follows... Figure 3 As shown in the figure, the confidence interval width expands significantly during the mid-heating stage as the yield of pyrolysis hydrocarbons increases rapidly, reflecting the direct impact of the uncertainty of pyrolysis kinetic parameters on the confidence level of the crack initiation pressure prediction.

[0144] In step S05, on day 10, the dynamic adjustment index is approximately 1.73, falling within the range of 1 to 2. The current time step is maintained, and Newton-Krylov iteration is employed. As heating continues, by day 24, the safety margin decreases to approximately 8.1 MPa, and the dynamic adjustment index is approximately 0.43, falling below 0.5. The system then reduces the time step to 0.25 times and triggers the warning process in step S06. The switching times of the dynamic adjustment index and iteration strategy at each point are shown in Table 2.

[0145] Table 2. Correspondence between dynamic adjustment index and iterative strategy at each critical moment.

[0146]

[0147] Step S06 is executed. On day 27, the safety margin drops to 4.8 MPa, below the warning threshold of 5 MPa, triggering a warning signal. The two-step approximation method calculates the current rate of change of fracturing pressure to be approximately... MPa / d, bottom hole pressure change rate approximately MPa / d, safety margin change rate approximately MPa / d; the estimated crack initiation time is approximately day 28.9. Based on this, the construction personnel will determine the supercritical... The injection pressure was reduced from 18 MPa to 14 MPa, and the heating rate was reduced from 5℃ / min to 3℃ / min. The recalculated crack initiation pressure rose in subsequent moments, and the safety margin was restored to above 7.2 MPa, allowing the construction to continue safely.

[0148] Compared with traditional methods for judging crack initiation, the advancement of this invention lies in the fact that traditional methods rely on a single elastic mechanics formula or empirical relationship, lack a quantitative description of the dynamic evolution of the hydrocarbon generation and pressurization process of kerogen pyrolysis over time, and cannot characterize supercritical processes. The continuous weakening effect caused by the dissolution of carbonate minerals leads to a systematic overestimation of the initiation pressure in the later stages of construction. This invention couples the pyrolysis kinetic equation with an elastic mechanics framework and embeds the spatial heterogeneity of pores and fractures and multi-field coupling constraints into the recursive prediction process using a physical constraint graph attention recursive network. This enables the initiation pressure prediction to respond in real time to the dynamic changes in pyrolysis and mineral weakening. The confidence interval provided by the Ritz adaptive approximation algorithm allows construction personnel to obtain quantitative prediction reliability information for the first time, rather than relying solely on single-point prediction values. The introduction of a dynamically adjusting exponential function automatically encrypts computational resources during dangerous construction phases and reasonably sparses them during safe phases, achieving an adaptive balance between accuracy and efficiency. This fundamentally compensates for the structural defects of traditional methods in multi-field dynamic coupling scenarios.

[0149] It should be noted that the variables involved in this invention are explained in detail in Tables 3, 4, and 5.

[0150] Table 3. Variable Explanation Table (Part 1)

[0151]

[0152] Table 4. Variable Explanation Table (Part Two)

[0153]

[0154] Table 5. Variable Explanation Table (Part 3)

[0155]

[0156] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for determining the initiation of cracks in low-maturity shale heated by supercritical carbon dioxide, characterized in that, Includes the following steps: By integrating conventional well logging, core experiments and formation testing data, the minimum horizontal principal stress, maximum horizontal principal stress, vertical stress, original formation pore pressure, formation temperature, total organic carbon content, carbonate mineral volume fraction, initial tensile strength, elastic modulus, Poisson's ratio, coefficient of thermal expansion, rock density, as well as supercritical carbon dioxide injection pressure, target heating temperature, heating time and heating rate parameters of medium and low maturity shale formations were obtained. Based on shale core samples of medium and low maturity, high temperature rock mechanics experiments, high temperature and high pressure closed hydrocarbon generation pressurization experiments, pyrolysis kinetics experiments, and mechanical comparison experiments after supercritical carbon dioxide immersion were carried out. The tensile strength temperature decay coefficient, pressurization conversion coefficient, discrete activation energy distribution model reaction kinetic parameters, strength weakening amplitude coefficient, and strength weakening rate coefficient were calibrated, and the effective porosity and pore throat connectivity coefficient were measured. The obtained parameters and calibration parameters are substituted into the multi-field coupled fracturing pressure formula to calculate the fracturing pressure at each time. Then, the physical constraint graph attention recursive network is input to output the corrected fracturing pressure and safety margin. The multi-field coupled fracturing pressure formula integrates the pressurization of kerogen pyrolysis hydrocarbon generation, the intensity reduction of supercritical carbon dioxide, the thermal expansion stress correction and the original geostress into the same analytical framework. Substituting the safety margin, the rate of change of the safety margin, and the uncertainty of the initiation pressure into the Ritz adaptive approximation algorithm, the spatial distribution of the initiation pressure and the confidence interval of the initiation pressure at each time point are output. Substitute the normalized values ​​of safety margin, safety margin change rate, and fracturing pressure uncertainty into the dynamic adjustment exponential function to calculate the dynamic adjustment exponential. Adjust the time step and iteration strategy of the physical constraint graph attention recursive network according to the interval to which the dynamic adjustment exponential belongs, and update the safety margin according to the difference between the bottom hole pressure and the fracturing pressure obtained from real-time monitoring. If the safety margin is lower than the warning threshold, the expected crack initiation time is obtained by linear extrapolation based on the safety margin and the rate of change of the safety margin using a two-step approximation method. The warning signal is then output and the construction parameters are adjusted.

2. The method according to claim 1, characterized in that, The minimum horizontal principal stress, maximum horizontal principal stress, and vertical stress are obtained through a combination of formation testing and well logging interpretation; the original formation pore pressure is obtained through formation pressure testing; and the formation temperature is obtained through thermometry logging.

3. The method according to claim 2, characterized in that, The calibration of the tensile strength temperature decay coefficient is specifically carried out by conducting Brazilian splitting tests on the same batch of core samples at 25℃, 100℃, 200℃, 300℃ and 400℃ respectively, and obtaining the temperature and rock tensile strength data by least squares fitting using an exponential decay model, with no less than 3 samples for each temperature point.

4. The method according to claim 3, characterized in that, The calibration of the reaction kinetic parameters of the discrete activation energy distribution model is specifically performed by thermogravimetric-mass spectrometry (TGA) experiments under multiple heating rates. The parameters are obtained by fitting the discrete activation energy distribution model with a Bayesian inference framework, which uses Markov chain Monte Carlo sampling to globally estimate the high-dimensional parameter space.

5. The method according to claim 4, characterized in that, The calibration of the strength weakening amplitude coefficient and the strength weakening rate coefficient is specifically carried out by immersing the core sample in a supercritical carbon dioxide environment for 1 day, 3 days, 7 days, 14 days and 30 days. After taking it out, the uniaxial compressive strength and rock tensile strength are measured, and the rock strength-time data are obtained by fitting the exponential decay model.

6. The method according to claim 5, characterized in that, In the multi-field coupled crack initiation pressure formula, the minimum and maximum horizontal principal stresses after temperature correction are calculated using the elastic constitutive relationship through the coefficient of thermal expansion, elastic modulus, and temperature change; the tensile strength of the rock after temperature correction is calculated using the exponential decay function through the temperature decay coefficient of tensile strength; and the pyrolysis hydrocarbon generation pressurization is calculated using the pressurization conversion coefficient, rock density, hydrocarbon yield, total organic carbon content, and effective porosity.

7. The method according to claim 6, characterized in that, The dynamic adjustment exponential function is a weighted sum of three data points: the normalized value of the safety margin, the normalized value of the rate of change of the safety margin, and the normalized value of the uncertainty of the crack initiation pressure.

8. The method according to claim 7, characterized in that, The physical constraint graph attention recursive network abstracts the shale pore and fracture network into a dynamic graph structure. The graph attention layer adopts a multi-head attention mechanism, introducing physical residual jump connections within each gated loop unit time step. The numerical residual of the pyrolysis kinetic equation and the conservation residual of the continuous equation are concatenated into an additional vector and superimposed on the candidate hidden state, driving the internal iterative correction of the network until the residual is lower than the physical residual convergence threshold.

9. The method according to claim 8, characterized in that, The training dataset of the physical constraint graph attention recursive network covers a sample set of parameter combinations under the conditions of total organic carbon content of 1% to 12%, carbonate mineral volume fraction of 5% to 40%, heating rate of 1℃ / min to 20℃ / min, and supercritical carbon dioxide injection pressure of 7MPa to 30MPa. It is divided into training set, validation set and test set in a ratio of 8:1:

1.

10. The method according to claim 9, characterized in that, The Ritz adaptive approximation algorithm is based on the Ritz variational principle of elasticity. It selects an adaptive basis function family that includes exponential thermal source term functions and Weibull distribution mineral weakening functions. Through the Rayleigh-Ritz process, it transforms the infinite-dimensional functional optimization into a finite-dimensional linear equation system. During the iteration process, it adaptively densifies or merges the basis function arrangement according to the residual distribution of the current solution.

Images