A grouting pressure and flow integrated optimization method based on diffusion inversion and gas grading

Abstract

The application provides a grouting pressure and flow integrated optimization method based on diffusion inversion and gas grading, and relates to the technical field of grouting pressure and flow optimization. It comprises the following steps: constructing a multi-field coupling prior model and, based on the final multi-field coupling prior model, collecting on-line field data and performing data preprocessing to obtain a time series data set; taking the final multi-field coupling prior model as a forward modeling framework and taking the time series data set as a constraint, constructing a joint objective function and determining a parameter posterior distribution; based on the parameter posterior distribution and gas state information, dynamically evaluating risks and performing multi-objective optimization, and outputting a joint control curve and a corresponding execution plan. The application solves the problem that, in the prior art, real-time closed-loop control cannot be realized under uncertain conditions to achieve parameter distinguishability, dynamic gas grading linkage and opportunity constraint multi-objective optimization.

Description

Technical Field

[0001] This invention relates to the field of grouting pressure and flow rate optimization technology, and in particular to an integrated optimization method for grouting pressure and flow rate based on diffusion inversion and gas classification. Background Technology

[0002] Current technologies generally employ constant pressure or constant flow grouting, using empirical interpretation of pressure-flow curves to infer the diffusion radius and grout absorption capacity of the surrounding rock. While mechanistic studies have established seepage-rheology coupling and equivalent permeability evolution models, their characterization of time-varying and spatial heterogeneity in the field remains insufficient. Regarding parameter inversion and data utilization, methods combining numerical forward modeling for parameter identification and preliminary data assimilation have emerged, capable of recovering key parameters such as equivalent permeability and yield stress to some extent. However, most of these methods remain offline, single-scale, and have limited sensing channels. Gas control has established regional classification and outburst prevention systems at the regulatory level; however, grouting control parameters are mostly set based on static levels or regional averages, lacking online classification and differentiated control coupled with the grouting diffusion state. In terms of optimized control, existing research often focuses on single objectives (such as injection volume or pressure threshold), with weak synergistic coordination of multiple objectives and constraints, including reinforcement radius uniformity, plugging efficiency, gas extraction gain, and safety probability.

[0003] To address the aforementioned limitations, the industry is evolving towards dynamic closed-loop and multi-scale integration: relying on multi-channel on-site sensing of pressure, flow rate, acoustic emission / microvibration, and gas pore pressure and concentration, continuous linkage between inversion and optimization is achieved through rolling time windows; at the model level, the transition is from single-scale to multi-scale assimilation of mechanism and data, using analytical / numerical forward modeling as the framework, and data-driven correction to improve robustness and adaptability; in boundary control, the shift is from fixed operating conditions to programmable excitation, actively shaping the diffusion front and improving parameter sensitivity and identifiability through methods such as pulsation, segmented step, staged pressure holding and depressurization; in gas management, the upgrade is from static zoning to dynamic hierarchical maps that integrate desorption dynamics and structural evolution, achieving real-time constraints coupled with grouting diffusion; in the optimization paradigm, the shift is from deterministic design to opportunity constraints and confidence-driven decision-making that consider the propagation of uncertainty, improving safety, interpretability, and engineering practicality.

[0004] However, several key bottlenecks still need to be overcome: First, the lack of tight coupling between inversion and optimization makes it difficult for control strategies to accurately match the actual diffusion front and gas state; second, the sensitivity of grouting observations under fixed boundaries to key parameters is limited, easily leading to slow inversion convergence, ambiguous results, and insufficient identifiability; third, gas classification and risk assessment have failed to achieve online linkage with the evolution of permeability and gas pressure field; fourth, there is a lack of an optimization framework that systematically weighs and updates objectives and constraints such as reinforcement quality, plugging efficiency, extraction gain, construction period and energy consumption, and safety probability; fifth, the decision-making process generally does not explicitly incorporate uncertainty propagation and confidence assessment, making it difficult to provide verifiable safety boundaries and reliable strategy switching logic. Summary of the Invention

[0005] To overcome the shortcomings of existing technologies, the purpose of this invention is to provide an integrated optimization method for grouting pressure and flow rate based on diffusion inversion and gas classification. This invention solves the problem that existing technologies cannot achieve real-time closed-loop control with parameter identification, dynamic gas classification linkage, and chance-constrained multi-objective optimization under uncertain conditions.

[0006] To achieve the above objectives, the present invention provides the following solution: An integrated optimization method for grouting pressure and flow rate based on diffusion inversion and gas classification includes: Collect multi-scale prior information and construct an initial multi-field coupled prior model; A grouting boundary excitation scheme is constructed, and the initial multi-field coupled prior model is parameterized with boundary excitation parameters to obtain the final multi-field coupled prior model. Based on the final multi-field coupling prior model, online field data is collected and data preprocessed to obtain a time-series dataset. Using the final multi-field coupled prior model as the forward modeling skeleton and the time series dataset as constraints, a joint objective function is constructed. Regularization and physical priors are introduced into the joint objective function to jointly estimate the spatiotemporally variable equivalent permeability, porosity, rheological parameters, interface leakage coefficient, and diffusion front geometry. The posterior distribution of the parameters is determined based on the results of the joint estimation. Based on the posterior distribution of parameters and gas state information, the risk is dynamically assessed and multi-objective optimization is performed, outputting the joint control curve of pressure, flow rate, and boundary excitation, along with the corresponding execution plan.

[0007] Preferably, the incentive scheme includes: Constant pressure, constant current, pulsation, segmented step, phased pressure holding and pressure release.

[0008] Preferably, the boundary excitation parameters include: Pulse frequency, duty cycle, step amplitude, pressure holding time, and pressure release rate.

[0009] Preferably, the field data includes: Pressure, flow rate, acoustic emission, displacement, temperature, gas pore pressure and concentration.

[0010] Preferably, the step of collecting multi-scale prior information and constructing an initial multi-field coupled prior model includes: Collect multi-source prior information, which includes: formation structure, mechanical and seepage parameters, initial gas pressure and content, fracture network and historical construction records; The prior information is processed by unifying units and coordinates, aligning time and space references, removing outliers and filling in missing information, and assigning prior uncertainty weights to obtain preprocessed multi-source prior information. The preprocessed multi-source prior information is mapped onto a unified geological model grid to obtain a spatially continuous initial parameter field; Based on the initial parameter field, establish the force, seepage, and gas coupling control equations and constitutive relations for grouting, surrounding rock, and gas, and set initial and boundary conditions to obtain the initial multi-field coupling prior model.

[0011] Preferably, the step of establishing the force, seepage, and gas coupling control equations and constitutive relations for grouting, surrounding rock, and gas based on the initial parameter field, and setting initial and boundary conditions to obtain an initial multi-field coupled prior model includes: The initial parameters are spatially mapped and interpolated on a unified grid to complete and smooth equivalent permeability, porosity, surrounding rock mechanical properties, fluid viscosity, rheological intensity and interface leakage capacity, thus obtaining a consistent parameter field. The constitutive relationships of grouting, surrounding rock, and gas are determined based on the aforementioned consistency parameter field; A set of governing equations is generated based on the constitutive relation; Initial and boundary conditions are set and the governing equations are coupled to establish an initial multi-field coupled prior model. The initial conditions include: the initial values ​​of the effective stress coupling coefficient, the initial values ​​of the liquid phase pressure, the initial values ​​of the liquid phase saturation, the initial values ​​of the equivalent permeability, the initial values ​​of the liquid phase equivalent viscosity, the initial values ​​of the gas phase saturation, and the initial values ​​of the gas phase pressure. Boundary conditions include: grouting orifice boundary excitation, surrounding rock mechanical boundary, seepage and gas boundary, interface and fracture boundary, compatibility and safety constraints, and monitoring and feedback interface. The grouting orifice boundary excitation is parameterized by boundary excitation parameters. The governing equations include: Mechanical control equations, seepage control equations, and gas control equations; The expression for the governing equations is as follows: ; The expression for the seepage control equation is: ; The expression for the gas control equation is: ; The expression for the initial multi-field coupled prior model is: ; in, For the initial multi-field coupled prior model, It is a displacement vector. The initial field for the displacement vector; The gradient is symmetric for the displacement; To follow the curing indication status The anisotropic stiffness tensor of evolution, The initial values ​​for the anisotropic stiffness tensor; The effective stress coupling coefficient affected by curing. The initial value for the effective stress coupling coefficient; and These are the initial values ​​of the liquid phase pressure and the liquid phase pressure, respectively. It is a unit second-order tensor; Body density; Porosity as it evolves with curing, and This represents the initial value for porosity; and This represents the liquid phase saturation and its initial value. The density of the liquid phase; This refers to the Darcy velocity in the liquid phase. For liquid phase source sinks; The equivalent permeability is modulated by both curing and anisotropy. This is the initial value for the equivalent permeability; and These are the relative permeability functions of the liquid phase and the gas phase, respectively. The equivalent viscosity of the liquid phase as it evolves with solidification. This is the initial value of the liquid phase equivalent viscosity; The liquid phase pressure gradient is g; g is the acceleration due to gravity. and This represents the gas phase saturation and its initial value. This refers to the gas phase density. The equivalent density of the solid framework or adsorption matrix; It is a monotonic function of the adsorption-desorption equilibrium loading that depends on the gas phase pressure; This refers to the Darcy velocity in the gas phase. For gas phase source and sink terms; This refers to the viscosity of the gas phase. and These are the initial values ​​for the gas phase pressure and the gas phase pressure, respectively. For computational domain, The computational domain boundary is defined by the excitation of the grouting orifice boundary, the surrounding rock mechanics boundary, the seepage and gas boundary, the interface and fracture boundary, the compatibility and safety constraints, and the monitoring and feedback interface. This is an initial constraint set that uniformly configures the excitation of mechanics, the seepage and gas boundaries, and the grouting orifice boundaries. The effective stress coupling coefficient is determined by the formation structure, mechanical and seepage parameters, fracture network, and historical construction records, while the liquid phase saturation is determined by the initial gas pressure and content.

[0012] Preferably, a grouting boundary excitation scheme is constructed, and the initial multi-field coupled prior model is parameterized with boundary excitation parameters to obtain the final multi-field coupled prior model, including: On the boundary of the initial multi-field coupled a priori model, the boundaries of grouting orifice, surrounding rock mechanics, seepage outer edge, gas outer edge, and fracture and interface are defined, and their geometric range and normal direction are marked to form a parameterizable set of boundaries. Based on the parameterizable boundary set, programmable basic elements integrating pressure and flow are defined for the grouting orifice boundary, standard basic elements of constant pressure, constant flow, and no flux are defined for the outer edge boundary, and displacement and stress basic elements are defined for the mechanical boundary. The parameterized interface of time history, amplitude and switching strategy is uniformly specified to obtain the boundary element library. Based on the boundary element library, the time history of the orifice pressure and flow rate is constructed and applied to obtain the grouting orifice time history excitation pair; Based on the time history excitation pair of the grouting orifice coupled with the fracture and matrix as a consistent expression of normal flux and pressure difference, the interface flux constraint set is obtained; Based on the interface flux constraint set, compatibility checks and scheduling are performed on the mechanical, seepage and gas boundaries to generate an executable boundary time history plan. The relationship between grouting orifice excitation and interface flux is embedded into the governing equations according to the boundary type, and the consistency of the initial conditions and parameter fields is checked to obtain a closed solution body containing time history boundaries. Based on the closed solution body containing the time history boundary, the operator form transformation of the boundary excitation parameterization is performed to obtain the final multi-field coupled prior model; Among them, the time history excitation pair of the grouting orifice includes: the time history of the grouting orifice pressure and the time history of the grouting orifice flow. The interface flux constraint set includes: the normal leakage flux density of the surrounding rock interface and the interface conductivity coefficient; The expressions for the grouting orifice time-history excitation pair, the interface flux constraint set, and the final multi-field coupled prior model are as follows: ; ; ; in, The time history of the grouting orifice pressure; The time history of the grouting orifice flow rate; and The baseline pressure and baseline flow rate for grouting; The pressure and increase in traffic during the phased leap; The pressure and flow rate amplitudes of the pulsating segment; The pulsation frequency; This represents the phase difference between pressure and flow rate. This represents the start time of the corresponding stage; The duration of the phase; A unit rectangular window function is used to define the phase window; denoted as the normal leakage flux density at the surrounding rock interface; The interface conductivity coefficient; The equivalent pore pressure in the interface neighborhood; This is the uncured flux base value; The directional conductivity coefficient in the fully cured state; The unit vector is the interface normal. It is a boundary operator with time and parameter vector as independent variables; For the boundary excitation parameter vector, This is the final multi-field coupled prior model after parameterization of boundary excitations; This is a computable solution body after embedding time history boundaries and completing consistency verification; A set for boundary classification and role labeling; For boundary primitive library; This is the coordinated time history set of grouting orifice pressure and flow rate; For the set of interface flux constraints; For multi-boundary compatible scheduling plans.

[0013] Preferably, based on the multi-field coupled prior model, online acquisition of field data and data preprocessing are performed to obtain a time-series dataset, including: Under a unified grid and geometric framework, the state variables and boundary positions of the pressure, flow rate, displacement, gas concentration and pore pressure monitoring channels are mapped one by one with the multi-field coupled prior model to obtain a structured online observation layout; Based on a structured online observation layout, raw data streams are continuously collected from each monitoring channel according to a predetermined sampling period and aligned and deduplicated using a unified time reference to obtain time-synchronized raw time-series records. The original time series records are subjected to quality inspection and saturation, jump, disconnection and drift anomalies are removed to obtain clean time series data; The cleaning time series data is subjected to zero-phase finite bandwidth filtering according to the physical bandwidth of the monitoring indicators and power frequency and mechanical harmonic interference is suppressed to obtain pre-filtered time series data; The pre-filtered time series data of each channel are resampled at a uniform target sampling rate, and short-time missing measurements are filled by unbiased interpolation based on local trends to obtain aligned time series data. The preprocessed synchronization segments from each monitoring channel are merged into a structured time-series dataset in chronological order.

[0014] Preferably, the final multi-field coupled prior model is used as the forward modeling skeleton, and the time-series dataset is used as a constraint to construct a joint objective function. Regularization and physical priors are introduced into the joint objective function to jointly estimate the spatiotemporally variable equivalent permeability, porosity, rheological parameters, interface leakage coefficient, and diffusion front geometry. The posterior distribution of the parameters is determined based on the joint estimation results, including: Under the framework of multi-field coupled forward modeling, an extended state vector containing field states and parameters to be estimated is constructed and time discretized and grid aligned to support the recursive solution of filtering, so as to obtain the extended state definition for assimilation. The parameters to be estimated include spatiotemporally variable equivalent permeability, porosity, rheological parameters, interface leakage coefficient and diffusion front geometry. Based on the time series dataset and monitoring layout, observation operators are defined to map the model state to the observation space, resulting in an observation vector flow and observation operators that match the extended state. The observation vector flow, observation operator, and spatial smoothing regularization term are used to construct a joint objective function and embed physical priors to estimate the parameters to be estimated, thereby obtaining the posterior estimates of the state and parameters at the current time. Based on the posterior estimates of the current state and parameters, the forward modeling skeleton is used to advance the extended state and covariance, and tangentization is used to process nonlinearities to obtain the prior distribution for the next time step. The prior distribution is corrected by observation assimilation and the state and parameters are jointly updated to obtain the posterior distribution at the next time step, and the updated posterior mean and covariance of the state and parameters are obtained. After each update, perform regularized consistent projection and physical domain correction to satisfy nonnegativity, upper bound, anisotropic symmetry and geometric consistency with the diffusion front, and obtain the posterior distribution of parameters; The extended state definition expression, the joint objective function expression, and the expression for the posterior parameter distribution are respectively: ; ; ; in, A timeline for boundaries and external stimuli; This is a preprocessed structured time-series dataset; For the first The extended state vector contains the field state. With the parameters to be estimated It is a nonlinear state transition operator; External stimulus and boundary input; The process noise has a covariance of For observation vectors; For observation operators; To observe the noise, the covariance is... The joint objective function; This is the parameter vector from the previous time step; For spatial difference operators; For spatial smooth covariance; The physical-consistent residual vector; For physical residual covariance; and These are the predicted and updated state estimates, respectively; and These are the predicted and updated covariances, respectively. The tangent matrix or set statistics are used to replace the state transition; Tangentization matrix or set statistics substitution for observation operators; Kalman gain; It is the identity matrix; and Let be the mean and covariance of the posterior parameters; and Let be the mean and covariance of the posterior state.

[0015] The present invention discloses the following technical effects: This invention provides an integrated optimization method for grouting pressure and flow rate based on diffusion inversion and gas classification, comprising: collecting multi-scale prior information and constructing an initial multi-field coupled prior model; constructing a grouting boundary excitation scheme, and parameterizing the initial multi-field coupled prior model with boundary excitation parameters to obtain a final multi-field coupled prior model; based on the final multi-field coupled prior model, collecting field data online and performing data preprocessing to obtain a time-series dataset; using the final multi-field coupled prior model as the forward modeling skeleton and the time-series dataset as constraints to construct a joint objective function, and introducing regularization and physical priors into the joint objective function to jointly estimate the spatiotemporally variable equivalent permeability, porosity, rheological parameters, interface leakage coefficient, and diffusion front geometry, and determining the posterior distribution of parameters based on the joint estimation results; dynamically assessing risks and performing multi-objective optimization based on the posterior distribution of parameters and gas state information, and outputting the joint control curve of pressure, flow rate, and boundary excitation and the corresponding execution plan. This invention addresses the problems in the background technology, such as parameter unidentification caused by constant pressure and constant flow, lack of robustness in point estimation parameter tuning, and difficulty in handling single-field model bias and field non-stationarity, through an online closed loop of "programmable boundary excitation - diffusion inversion - uncertainty-driven optimization". It proposes an integrated solution: on the one hand, it improves the identifiability of equivalent permeability, porosity, rheological parameters, interface leakage coefficient and diffusion front geometry by using boundary excitations such as pulsation, step, and pressure holding / releasing, significantly reducing posterior variance and multiple solutions; on the other hand, it performs rolling assimilation of multi-source time series data under the unified force-permeability-gas forward model framework, maintaining real-time updates of parameters and states, and shortening the test injection and tuning cycle; furthermore, it introduces the posterior distribution of parameters and gas state information into multi-objective optimization and chance-constrained control, so as to simultaneously meet safety and efficiency objectives under pre-set confidence, and finally outputs the joint control curve and execution plan of pressure-flow-boundary excitation, improving robustness, efficiency and reproducibility. Attached Figure Description

[0016] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0017] Figure 1 A flowchart of an integrated optimization method for grouting pressure and flow rate based on diffusion inversion and gas classification is provided for an embodiment of the present invention. Figure 2 A flowchart for constructing an initial multi-field coupling prior model provided in this embodiment of the invention; Figure 3The flowchart for constructing the final multi-field coupling prior model provided in the embodiments of the present invention is shown. Detailed Implementation

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

[0019] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0020] like Figure 1 As shown, this invention provides an integrated optimization method for grouting pressure and flow rate based on diffusion inversion and gas classification, comprising: Step 100: Collect multi-scale prior information and construct an initial multi-field coupled prior model; Step 200: Construct a grouting boundary excitation scheme, and parameterize the initial multi-field coupled prior model with boundary excitation parameters to obtain the final multi-field coupled prior model; Step 300: Based on the final multi-field coupled prior model, collect field data online and perform data preprocessing to obtain a time-series dataset; Step 400: Using the final multi-field coupled prior model as the forward modeling skeleton and the time series dataset as constraints, construct a joint objective function, and introduce regularization and physical priors into the joint objective function to jointly estimate the spatiotemporally variable equivalent permeability, porosity, rheological parameters, interface leakage coefficient and diffusion front geometry, and determine the posterior distribution of parameters based on the joint estimation results; Step 500: Based on the posterior distribution of parameters and gas state information, dynamically assess the risk and perform multi-objective optimization, outputting the joint control curve of pressure, flow rate, and boundary excitation, and the corresponding execution plan.

[0021] Furthermore, such as Figure 2As shown, the implementation process of step 100 is as follows: In this embodiment, the multi-source prior information is first collected in a structured manner and subjected to constraint verification to ensure data consistency and traceability. For the stratigraphic structure and fracture network, a discrete set of elements of stratigraphic position, fault, fracture orientation and connectivity is formed by seismic interpretation, borehole logging, well logging and geological cataloging, and geometric legality check and topological closure are performed using a unified coordinate benchmark and unit scale threshold. For mechanical and seepage parameters, indoor triaxial and permeability tests, field reinjection and drainage tests, and historical grouting comparison tests are collected to provide the measured values ​​and confidence intervals of elastic modulus, Poisson's ratio, yield strength, equivalent permeability, porosity and interfacial conductivity. For the initial gas pressure and content, gas extraction, pore pressure monitoring, gas content measurement and temperature field are integrated to form a spatial discrete observation set and mark the timestamp, sampling method and quality level. This embodiment applies coordinate integration, unit unification, and time reference alignment to all raw data. It also removes saturation, disconnection, and instrument drift anomalies using three types of constraints based on sensor range, gradient threshold, and steady-state drift detection. Local kriging or radial basis functions are used to fill in the missing data unbiasedly. At the same time, a priori uncertainty weight is assigned to each source data for subsequent weighted fusion and constraint propagation.

[0022] This embodiment implements physical consistency mapping and hierarchical fusion of prior information on a unified geological model grid to ensure the continuity of the parameter field and its compatibility with structural boundaries. For discrete fracture elements, this embodiment generates principal fields of anisotropic principal axes and equivalent permeability tensors through strike density and aperture statistics, and superimposes interface conductivity enhancement factors at fault and interlayer interfaces. For mechanical parameters, this embodiment sets constant or slowly varying fields according to lithological and structural domains, and balances porosity and equivalent permeability through volume constraints and empirical relationships. For the initial gas field, this embodiment uses pore pressure observation as a hard constraint and gas content as a soft constraint, and employs covariance kernel and stratigraphic zonal prior structural joint kriging to obtain the spatially continuous distribution of initial gas pressure and liquid phase pressure. Local smoothing and interface transition limiting are applied in the wellbore and interface neighborhood to form a consistent initial parameter field. To avoid loss of identifiability due to excessive smoothing, this embodiment applies high gradient edge-preserving terms to fracture zones and fault neighborhoods to retain key heterogeneous features.

[0023] This embodiment establishes the force, seepage, and gas coupling control equations and constitutive relationships for grouting, surrounding rock, and gas based on a consistent initial parameter field, and clarifies the initial and boundary conditions to form a computable initial multi-field coupled a priori model. The mechanical side adopts effective stress and anisotropic linear elastic or piecewise elastoplastic constitutive models; the seepage side adopts two-phase Darcy and introduces relative permeability functions and interface conductivity constraints; and the gas side adopts isothermal gas state and kinetic closure of adsorption loading. This embodiment provides the initial spatial distribution of displacement, liquid phase pressure, gas phase pressure, liquid and gas phase saturation, porosity, equivalent permeability, liquid phase viscosity, effective stress coupling coefficient, and adsorption loading at the start of the calculation. Compatible boundary conditions and time-history constraints are configured according to the wellhead boundary, outer edge boundary, interface and fracture boundary, mechanical support boundary, and monitoring feedback interface to ensure conservation, symmetry, and stability requirements. Mature technologies such as data format conversion, conventional coordinate projection, standard numerical discretization, and solution processes are not elaborated upon.

[0024] Specifically, this embodiment performs spatial mapping and interpolation of initial parameters on a unified grid. First, a geological model grid is established using the same coordinate reference and unit topology. Discrete measurement values ​​such as permeability, porosity, elastic modulus, Poisson's ratio, yield strength, liquid viscosity, and interfacial conductivity obtained from field and laboratory measurements are interpolated intradomain and fused across domains according to lithological zones and structural boundaries. Intradomain interpolation uses a covariance kernel based on correlation length and anisotropic principal axes. Data-free areas are filled using soft-constrained priors consistent with stratigraphic zones. Cross-domain fusion sets boundary-preserving terms at faults and interlayer interfaces to maintain the physical rationality of parameter transitions. To avoid excessive smoothing leading to reduced identifiability, gradient limiting and boundary fidelity are applied to fracture-rich zones to preserve key heterogeneity.

[0025] Based on the above-mentioned consistent parameter field, the constitutive relationship between grouting, surrounding rock, and gas is determined: the surrounding rock adopts anisotropic linear elastic or segmented elastoplastic constitutive model, and the parameters are given by the inversion values ​​of indoor triaxial and field reinjection tests; the seepage side adopts two-phase Darcy flow, and the relative permeability functions of the liquid phase and gas phase are fitted according to the saturation curve and extrapolated by the empirical curves of adjacent lithologies in the unmeasurable area; the interface conductivity is characterized by a monotonic function that changes with the solidification indication state, and the conductivity baseline value is taken when it is not solidified, and it converges to the low conductivity limit along the anisotropic principal axis after complete solidification. The governing equations are generated sequentially according to mechanics, seepage, and gas, and initial and boundary conditions are configured at the start of the calculation. The initial conditions include the effective stress coupling coefficient, liquid phase pressure, gas phase pressure, liquid and gas phase saturation, equivalent permeability, and spatial distribution of liquid phase viscosity, which are derived from pre-processed multi-source prior information and the matching state equations and empirical relationships. The boundary conditions include the pressure and flow excitation of the grouting orifice over time, the displacement or stress constraint of the surrounding rock, the constant pressure or no-flux condition of the seepage and gas outer edge, the normal leakage constraint of the fracture and interface, and the observation mapping of the monitoring and feedback interface. The grouting orifice excitation is obtained through parameterized configuration, and its parameters include baseline amplitude, step increment, pulsation amplitude, frequency, phase, stage start and end time, and duration, forming a programmable time history.

[0026] For ease of understanding, the consistency parameter field refers to the set of spatially continuous parameters obtained after interpolation, completion, and smoothing under a unified grid, unified coordinates, and physical constraints. It is used to ensure that each governing equation is numerically closed and consistent across fields. The solidification indicator state refers to the degree of material solidification accumulated at pores and interfaces over time and during construction. Its value ranges from unsolidified to fully solidified, determining the extent of stiffness enhancement, viscosity increase, and interface conductivity reduction. Interface conductivity refers to the equivalent passage capacity of fluid flux along the normal and tangential directions of fractures or phase interfaces. Its value range is calibrated based on the on-site seepage pressure response and historical grouting records. In the numerical implementation, the parameter sources, values, and functions of each constituent relationship and control term are explained as follows: The effective stress coupling coefficient is jointly given by core tests, experience with similar formations, and historical construction records, and is used to convert changes in pore pressure into adjustments to the effective stress of the surrounding rock; the initial values ​​of pressure and saturation of the liquid and gas phases are estimated by pore compaction measurements, gas content, and temperature fields, and are generally obtained through a joint inversion of near-well hard constraints and far-field soft constraints; the initial values ​​of equivalent permeability and porosity are obtained based on indoor permeability tests, field reinjection tests, and fracture statistics, and are decomposed into principal values ​​according to the anisotropic principal axis. The direction set; the initial value of the liquid phase equivalent viscosity is determined by the cement slurry or chemical slurry ratio and temperature, and increases monotonically with the solidification indication state; the relative permeability function is obtained by fitting the experimental saturation curve, and is used to allocate the relative passage capacity of the two phases in the pores; the interface leakage parameters are given by the on-site calibration of pressure drop and flux, and decrease with the solidification indication state to reflect the sealing effect; the volume density and gravity direction are determined according to the engineering site agreement; the monitoring and feedback interface maps the pressure, flow, displacement and gas concentration sensing channels into observation points and observation time series, which are used for subsequent verification of the consistency between the initial conditions and the boundary time history.

[0027] Furthermore, such as Figure 3 As shown, the implementation process of step 200 is as follows: This embodiment achieves boundary excitation parameterization and obtains the final multi-field coupled prior model through four steps: boundary partitioning, primitiveization, temporal history synthesis, and consistency verification. First, the boundaries are classified and labeled within a unified geometry and mesh framework to form a set of parameterizable boundaries, specifically including grouting orifice boundaries, surrounding rock mechanical boundaries, seepage and gas outer edge boundaries, and fracture and interface boundaries. The geometric range and normal direction of each type of boundary are recorded. Secondly, a boundary element library is constructed. The grouting orifice boundary is configured with integrated programmable elements for pressure and flow. The outer boundary provides three types of standard elements: constant pressure, constant flow, and no flux. The mechanical boundary provides displacement and stress elements. All elements uniformly expose parameterized interfaces for time history, amplitude, and switching strategies. Furthermore, based on the boundary element library, the coordinated time history of grouting orifice pressure and flow rate is synthesized, specifically formed by splicing a baseline segment, a segmented step segment, and a pulsating segment. The baseline segment is used to set the grouting initiation level, the segmented step segment is used to trigger an identifiable diffusion response, and the pulsating segment is used to improve the identification of equivalent permeability and interface leakage. The custom term "grouting orifice time history excitation pair" refers to the paired input consisting of pressure time history and flow rate time history. At the same time, the fracture and matrix interface is uniformly described as a relationship where the normal flux is driven by pressure difference and modulated by the interface conductivity coefficient. The custom term "interface conductivity coefficient" refers to a parameter reflecting the interface throughput capacity, which can monotonically transition from the higher throughput capacity of uncured to the lower throughput capacity of fully cured, and allows the directionality to be reflected by the angle between the interface normal and the anisotropic principal axis. Subsequently, compatibility checks and scheduling are carried out to ensure that the coordinated time history of orifice pressure and flow does not conflict with the constant pressure or no flow at the outer edge, mechanical displacement or stress conditions, and interface flux constraints within the same time history. If necessary, contradictions are eliminated by redistributing the duty cycle and stage duration to form an executable boundary time history plan. Finally, the relationship between the grouting orifice excitation and the interface flux is written into the boundary condition terms of the governing equations according to the boundary type. The consistency of the initial conditions and parameter field is then checked. This check includes matching the initial pressure and saturation with the constant pressure or no flux at the outer edge, the continuity and upper limit of the orifice stage switching, and the consistency and monotonicity of the interface flow parameters and solidification indicators. After verification, a closed solution body containing time-history boundaries is obtained. Based on this, an operator-based transformation of the excitation parameterization is performed. A custom boundary operator is defined, taking time and boundary excitation parameters as independent variables, and mapping the pressure and flow time histories, interface flux constraints, and boundary classifications into a set of boundary inputs that can directly drive the solution of the governing equations. The parameters of the object include baseline pressure and baseline flow rate, segmented step amplitude, pulsation amplitude, pulsation frequency, phase difference, stage start time, stage duration, duty cycle, pressure holding time, and pressure release rate. Among them, the baseline parameters can be determined by the equipment calibration and formation allowable pressure before the test injection. The step and pulsation amplitudes are set based on identifiability analysis and safety upper limits. The frequency is selected based on the formation relaxation time constant to avoid being too fast or too slow. The stage time is determined by the construction plan and available monitoring time window. The duty cycle is used to control the effective power of the pulsation. The pressure holding and pressure release parameters correspond to the grouting quality and safe release requirements. The initial value of the interface flow parameters is calibrated by historical grouting and pressure drop tests and is monotonically updated with the solidification process. Through the above implementation process, the boundary excitation schemes of constant voltage, constant current, pulsation, segmented step, and staged holding or releasing of pressure can all be combined and called as programmable primitives, and finally output a reproducible boundary time history plan and a parameterized multi-field coupled prior model.

[0028] More specifically, the "unit rectangular window function" in this embodiment is used to precisely define the effective interval of a certain excitation stage on the time axis. This means that it takes a value of one during the time period covered by the preset stage start time and duration, and a value of zero at other times, thereby achieving segmented start-stop switching control of the pressure and flow time histories. Its parameters are given by the construction plan and equipment response capabilities, including the stage start time and stage duration, and can be combined with duty cycle and phase difference to construct segmented pulsations, pressure holding, or step segments. In implementation, the target segment shape of pressure or flow is multiplied by the window function and then summed according to adjacent time segments to form a coordinated time histories. Furthermore, aligning the window boundary with the solution time step avoids amplitude leakage and phase shift. This window function acts simultaneously on both pressure and flow channels, ensuring that the "grouting orifice time histories excitation pair" are coordinated within the same time window, and together with compatibility and safety constraints, ensures the continuity and safety of inter-segment switching.

[0029] Furthermore, the implementation process of step 300 is as follows: This embodiment first completes the structured mapping and calibration of the online observation layout under a unified grid and geometric framework: based on the state variables and boundary classification in the multi-field coupled prior model, pressure and flow sensors are positioned to the corresponding units or face centers of the grouting orifice boundary, displacement gauges and strain gauges are mapped to the nodes or sections of the surrounding rock mechanical boundary, and gas pore pressure gauges and gas concentration probes are mapped to the volume unit centers of the seepage and gas outer edge boundaries or the borehole-tunnel neighborhood, respectively, and their normal direction, sampling surface or control volume geometric scale is recorded to determine the numerical observation operator; for acoustic emission channels, this embodiment arranges and records the sensitive volume and frequency band range in the neighborhood of fracture and interface boundaries for subsequent event-level matching with the micro-fracture activity index in the model; the above mapping generates a unique index from the measuring point to the model unit, forming an observation metadata list containing channel type, spatial coordinates, normal direction, controlled surface / volume scale, corresponding state variables and boundary roles, which serves as a fixed reference for data acquisition and assimilation.

[0030] After obtaining the structured layout, this embodiment performs online acquisition and time series preprocessing: each monitoring channel is triggered according to a predetermined sampling period and a unified time reference. The sampling side adopts a time stamp alignment and deduplication strategy, prioritizing the retention of the earliest timestamp within the trigger window that has passed verification, forming a time-synchronized original time series. For this original time series, this embodiment implements a quality inspection process, setting saturation upper and lower limits automatically derived according to the device range. The jump criterion adopts the joint identification of adjacent sample difference and adaptive median absolute deviation threshold. Disconnection is identified by heartbeat absence and long silent window. Drift is identified by sliding window linear trend and low-frequency energy ratio. All abnormal points are marked with their source and removed or downweighted to obtain clean time series data. Subsequently, zero-phase finite bandwidth filtering is implemented according to the effective bandwidth of each physical quantity. The pressure and flow channels set notches at the power frequency and its harmonics and retain the low-frequency to mid-frequency components required for diffusion and pulsation identification. The displacement channel suppresses high-frequency mechanical vibration. The gas pore pressure and concentration channels suppress pumping cycle harmonics. The acoustic emission channel retains its working bandwidth and strongly suppresses DC and power frequency components, obtaining pre-filtered time series data for each channel.

[0031] This embodiment concludes with resampling, completion, and structured merging to form a time-series dataset suitable for assimilation: Pre-filtered data from each channel are resampled at a uniform target sampling rate, selected based on the highest physical bandwidth and Nyquist condition while also considering computational load. For short-term missing data, this embodiment employs unbiased interpolation based on local trends. For single-peak or step-adjustment neighborhoods, piecewise linear and local splines are combined with extrapolation suppression. Phase and amplitude are preserved without distortion for pulsating segments. For segments with continuous missing data exceeding the threshold, only missing data markers and uncertainty weighting are applied without interpolation. Subsequently, all channels are segmented into synchronous segments on the aligned time axis, edge samples at discontinuous points are removed, and the data is merged in chronological order into a structured time-series dataset. This dataset includes observation metadata, quality flags, filtering and resampling parameter records, and uncertainty weights, used for interfacing with the observation operators of multi-field coupled prior models to achieve subsequent rolling assimilation and boundary excitation verification.

[0032] Furthermore, the implementation process of step 400 is as follows: This embodiment uses the final multi-field coupled prior model as the forward modeling framework to construct an extended state containing the field state and the parameters to be estimated. Time discretization and grid alignment are performed to support recursive solution and online assimilation. The extended state consists of a displacement field, a pore pressure field, a gas phase pressure field, a two-phase saturation field, and coupled spatiotemporal variable parameters. These parameters include the principal value and direction of the equivalent permeability, porosity, rheological parameters, interfacial leakage coefficient, and the geometric descriptor of the diffusion front. The diffusion front geometry refers to the location of the movable boundary or isosurface reflecting the influence range of grouting or seepage, which can be described by the discrete values ​​of the front arrival time field or the front indicator function on the grid. To ensure a one-to-one correspondence with the time-series dataset, this embodiment defines observation operators based on the observation layout, mapping model states and parameters to the observation space of each channel to form a time-synchronized observation vector flow. The pressure and flow channels are mapped to the state or flux at the orifice boundary and outer edge boundary; the displacement channel is mapped to the displacement or stress-strain mapping at the surrounding rock mechanical boundary; the gas pore pressure and concentration channel is mapped to the pressure and phase composition at the center of the volume element; and the acoustic emission channel is mapped to the event intensity and frequency band energy indices in the fracture zone neighborhood. These extended states and observation operators are bound to unique indices on a unified grid to ensure spatiotemporal alignment and repeatable computation, and the physical dimensions, numerically stable value range, and hard constraints of each type of variable are recorded for subsequent regularized consistent projection and physical domain correction.

[0033] This embodiment uses the observation vector flow, observation operator, and spatial smoothing regularization term to form a joint objective, and embeds physical priors to constrain the identifiability and physical rationality of the parameters to be estimated. The joint objective consists of three parts, corresponding to the observation residual penalty, spatial smoothing regularization, and physical prior stigmatization penalty, respectively. The observation residual is used to measure the difference between the model output and the preprocessed structured time-series dataset. The spatial smoothing regularization is used to suppress unnecessary high-frequency fluctuations and maintain physical continuity in the anisotropic principal axis direction. The physical prior stigmatization is used to constrain the evolution of diffusion front geometry and penetration anisotropy near the solidification process, fracture orientation, and interlayer interfaces. To ensure full disclosure, this embodiment explains the source and use of parameters involved in the joint objective as follows: the time history plan of the boundary and external excitation comes from the boundary excitation parameterization stage, providing time-driven input for the forward modeling skeleton; the structured time series dataset comes from the online acquisition and preprocessing stage, and is the target data of the observation residuals; the process noise and its covariance are used to quantify the model uncertainty of the forward modeling skeleton during off-walking advancement, and the baseline level can be given according to lithological zoning and appropriately adjusted upwards during the construction stage; the observation noise and its covariance are used to express the measurement error and preprocessing residual interference of each channel, and its diagonal elements can be set according to the equipment accuracy and bandwidth and updated over time; the spatial difference operator and spatial smoothing covariance are used to achieve anisotropic smoothing, and the principal axis direction is given by fracture and bedding statistics; the physical residual vector and its covariance are used to quantify the deviation from physical laws, such as the geometric consistency of the diffusion front and pressure isosurface, and the consistent monotonicity of interface leakage and solidification indicators; after minimizing the joint objective, the posterior estimate of the current state and parameters and its covariance are obtained. For custom terms, extended state refers to expressing the state field and the parameters to be estimated in the same vector to support joint updates; physical prior refers to the statistical constraint terms constructed from empirical laws and constraint relationships; and diffusion front geometry refers to the set of geometric parameters of the boundary or transition zone of the seepage or grouting influence range.

[0034] This embodiment employs a recursive framework of prediction and updating to advance the posterior distribution, ensuring consistency between nonlinear computability and physical constraints. In the prediction step, this embodiment advances the mean of the extended state using a forward modeling skeleton and approximates the covariance using ensemble statistics or tangent matrix approximation. Boundaries and external stimuli serve as time history inputs to drive state evolution, and process noise covariance is introduced to reflect unmodeled or discretized errors. In the update step, this embodiment uses observation vectors and observation operators to calculate gains, incorporates observation information into the prior distribution to obtain updated posterior mean and covariance, and employs localization and channel weighting when necessary to suppress long-range spurious correlations and reflect the reliability differences between different channels. After each update, this embodiment performs regularized consistent projection and physical domain correction, specifically including applying non-negativity and upper bound constraints to the equivalent permeability, porosity, interface leakage coefficient and rheological parameters, applying symmetry and principal axis order constraints to the anisotropic parameters, applying geometric monotonicity and topological connectivity constraints consistent with the pressure field and saturation field to the diffusion front geometry, and retaining the necessary high gradient in the neighborhood of the crack and interface to avoid recognition loss caused by excessive smoothing. To ensure feasibility, this embodiment illustrates the source and values ​​of key matrices and vectors with examples: tangentization or set statistical substitution can be constructed using unscented transformation or set perturbation; the initial covariance value is set according to similar engineering experience and adaptively adjusted during rolling assimilation; the Kalman gain is jointly determined by the prior covariance, the observation operator, and the observation noise covariance; the identity matrix is ​​used as a dimension-consistent identity mapping to stabilize numerical operations; the mean and covariance of the parameter posterior and state posterior are used as output results to generate the prior input for the next step; through the above recursion, the spatiotemporally variable equivalent permeability, porosity, rheological parameters, interface leakage coefficient, and diffusion front geometry parameter posterior distributions are obtained, achieving joint estimation and uncertainty quantification consistent with the time-series dataset.

[0035] More specifically, in this embodiment, "model state" refers to the set of all field quantities and their boundary derivatives obtained from the forward modeling skeleton on a unified grid and time step, which are available for observation mapping. Specifically, this includes the displacement field and the strain and stress that can be calculated from it on the rock mechanics side; the liquid phase pressure field, gas phase pressure field, liquid and gas phase saturation field, and the derived relative permeability and equivalent viscosity on the seepage and gas side; the normal flux, tangential flux, and pressure difference at the boundaries and interfaces that can be calculated from the internal domain state; and the spatiotemporal variable parameters that participate in coupling and need to be assimilated on the grid. The field values ​​on the grid, such as the equivalent permeability principal value and direction, porosity, rheological parameters, interface leakage coefficient, and geometric indicators of the diffusion front used to describe the range of influence, are used to extract and interpolate the components that are isomorphic to the physical quantities of each sensing channel in the above states to the sensor position or project them to the measurement surface. The unit and dimension conversion and direction projection are completed according to the channel response model, thereby generating an observation vector flow with the same dimension as the extended state at each time step and giving the corresponding observation operator, realizing a reproducible mapping from "internal model quantities" to "observable quantities".

[0036] Furthermore, the implementation process of step 500 is as follows: This embodiment first constructs a dynamic risk assessment framework based on the posterior distribution of parameters and gas state information, clarifying risk indicators, trigger thresholds, and calculation time windows. Risk indicators include the probability of surrounding rock instability, the probability of overpressure rupture, the probability of sudden inrush and escape, the risk of blockage caused by relative permeability imbalance, and the risk of interface leakage exceeding limits. All indicators use the mean and covariance of the posterior distribution of parameters as input, and are updated continuously based on the latest synchronous observations of current gas pressure, gas concentration, temperature, and pore pressure. This embodiment calculates the conditional risk curve for each indicator within the future rolling prediction time window, using Monte Carlo sampling or ensemble perturbation to generate multiple possible trajectories and statistically analyzing the frequency of exceeding the safety threshold, forming a time series of risk exceedance rate and remaining safety margin. Simultaneously, a compatibility constraint set is constructed, consisting of the equipment reachability range, boundary excitation upper limit, interface flow convergence conditions dependent on the solidification process, the target mortar or slurry viscosity window, the allowable orifice pressure drop rate, and the gas drainage linkage strategy. To ensure feasibility, this embodiment provides the source and update rules for each threshold. For example, the surrounding rock instability threshold is given by the rock strength parameter and the current effective stress inversion, the inrush threshold is determined by the gas pressure gradient and connectivity index, and the interface leakage threshold is given by the solidification degree indicator and compatibility test curve. It also allows for adaptive correction within a limited range based on monitoring feedback.

[0037] This embodiment constructs a multi-objective optimization problem based on risk assessment and provides specific executable definitions for decision variables, objective functions, and constraints. Decision variables are orifice pressure, orifice flow rate, and boundary excitation mode parameters at discrete-time grid points. Boundary excitation mode parameters include pulsation frequency, duty cycle, stage start and end times, piecewise step amplitude, pressure holding duration, and pressure release rate, allowing switching between constant pressure, constant flow, pulsation, piecewise step, and staged pressure holding or pressure release within a defined set. The objective function adopts a multi-objective form combining weighted and sequential tolerance, with the main objectives being to maximize injection efficiency and diffusion front coverage uniformity, minimize unit volume injection energy consumption and interface leakage loss, and minimize risk exceedance rate and uncertainty volume. A penalty term for parameter posterior entropy is introduced to encourage improved identifiability. Constraints include hard constraints from compatibility and safety constraint sets and logical constraints for operating condition switching. Logical constraints require mode changes to be completed within the allowed switching window, with amplitudes continuous before and after the switch or gradually changing within a specified slope. This embodiment adopts a two-stage solution strategy. First, a set of safe and feasible initial control curves are quickly generated within the narrowed feasible domain. Then, the initial solution is refined by rolling time-domain optimization. During the refinement process, risk-sensitive adaptive weighting is implemented on the posterior covariance of the parameters. That is, the excitation intensity is appropriately reduced or the holding time is shortened during periods of high uncertainty to reduce tail risk. At the same time, the effective power of the pulsating component is increased within the identifiable window to improve the parameter assimilation quality.

[0038] Finally, this embodiment translates the optimized joint control curve of pressure, flow rate, and boundary excitation into an executable plan, and provides execution, monitoring, and rollback strategies. The execution plan lists the target pressure, target flow rate, excitation mode and parameters, allowable deviation band, maximum slope, and triggering conditions for each stage in a timeline. Triggering conditions include monitoring quantity exceeding threshold, risk exceedance rate increase, equipment anomaly, and data quality degradation, etc. A priority queue is defined to handle concurrent triggers, ensuring that safety-related triggers have the highest priority. To ensure closed-loop operation, this embodiment configures a rolling update cycle. In each update cycle, the posterior distribution of parameters and risk indicators are refreshed using the latest structured time-series dataset. The control curve in the remaining time domain is re-optimized without changing the executed segments. If necessary, a rollback strategy is initiated. The rollback strategy prioritizes switching to a low-energy constant-pressure or constant-current conservative operating condition, increasing the depressurization rate or reducing the duty cycle. At the same time, observations are intensified at key locations on the interface and outer edge. During execution, all instructions, actual responses, deviations, and triggering events are recorded to form audit and retraining data. This data is used for subsequent recalibration of optimization weights and iterative updates of empirical thresholds, ensuring that the joint control curve and execution plan achieve an optimal reproducible trade-off between safety, efficiency, and identifiability.

[0039] Specifically, the gas state information in this embodiment is jointly obtained from online monitoring and model inference and is updated on a rolling basis under a unified time reference. This includes real-time measurements directly collected by gas pore pressure gauges and gas concentration probes deployed in orifice, tunnel, and surrounding rock units, followed by a synchronized sequence after noise reduction, notch filtering, and resampling as described above. Combined with short-term forward modeling driven by the current boundary excitation and posterior mean of physical properties using a multi-field coupled prior model, state interpolation and extrapolation for undeployed locations and short-term missing measurement intervals are obtained. Finally, the uncertainty metric provided by the joint assimilation results is used to weight and fuse the above measurements and predictions, forming a spatial overlay. This embodiment generates time-series snapshots of the gas pressure field, gas partial pressure field, and volume fraction field in the vicinity of the orifice, the fracture zone, and the air side of the roadway. It also calculates derived quantities including gas pressure differential, concentration gradient, ventilation flux, and connectivity indication. To ensure feasibility, this embodiment incorporates sensor calibration coefficients, installation depth and orientation, sampling period, quality indicators, and the response coefficients of the observation-model mapping into the metadata list. Furthermore, it stipulates that when data quality degrades or missing data exceeds a threshold, model predictions and trends at adjacent measuring points are replaced by a weighted average of uncertainty, ensuring the continuity and traceability of gas state information in time and space, and providing real-time input for risk assessment and multi-objective optimization.

[0040] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on the differences from other embodiments. The same or similar parts between the various embodiments can be referred to each other.

[0041] This document uses specific examples to illustrate the principles and implementation methods of the present invention. The descriptions of the above embodiments are only for the purpose of helping to understand the method and core ideas of the present invention. Furthermore, those skilled in the art will recognize that, based on the ideas of the present invention, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of the present invention.

Claims

1. A method for integrated optimization of grouting pressure and flow rate based on diffusion inversion and gas classification, characterized in that, include: Collect multi-scale prior information and construct an initial multi-field coupled prior model; A grouting boundary excitation scheme is constructed, and the initial multi-field coupled prior model is parameterized with boundary excitation parameters to obtain the final multi-field coupled prior model. Based on the final multi-field coupling prior model, online field data is collected and data preprocessed to obtain a time-series dataset. Using the final multi-field coupled prior model as the forward modeling skeleton and the time series dataset as constraints, a joint objective function is constructed. Regularization and physical priors are introduced into the joint objective function to jointly estimate the spatiotemporally variable equivalent permeability, porosity, rheological parameters, interface leakage coefficient, and diffusion front geometry. The posterior distribution of the parameters is determined based on the results of the joint estimation. Based on the posterior distribution of parameters and gas state information, the risk is dynamically assessed and multi-objective optimization is performed, outputting the joint control curve of pressure, flow rate, and boundary excitation, along with the corresponding execution plan.

2. The integrated optimization method for grouting pressure and flow rate based on diffusion inversion and gas classification as described in claim 1, characterized in that, The incentive scheme includes: Constant pressure, constant current, pulsation, segmented step, phased pressure holding and pressure release.

3. The integrated optimization method for grouting pressure and flow rate based on diffusion inversion and gas classification as described in claim 1, characterized in that, Boundary excitation parameters include: Pulse frequency, duty cycle, step amplitude, pressure holding time, and pressure release rate.

4. The integrated optimization method for grouting pressure and flow rate based on diffusion inversion and gas classification as described in claim 1, characterized in that, The on-site data includes: Pressure, flow rate, acoustic emission, displacement, temperature, gas pore pressure and concentration.

5. The integrated optimization method for grouting pressure and flow rate based on diffusion inversion and gas classification as described in claim 1, characterized in that, The process of collecting multi-scale prior information and constructing an initial multi-field coupled prior model includes: Collect multi-source prior information, which includes: formation structure, mechanical and seepage parameters, initial gas pressure and content, fracture network and historical construction records; The prior information is processed by unifying units and coordinates, aligning time and space references, removing outliers and filling in missing information, and assigning prior uncertainty weights to obtain preprocessed multi-source prior information. The preprocessed multi-source prior information is mapped onto a unified geological model grid to obtain a spatially continuous initial parameter field; Based on the initial parameter field, establish the force, seepage, and gas coupling control equations and constitutive relations for grouting, surrounding rock, and gas, and set initial and boundary conditions to obtain the initial multi-field coupling prior model.

6. The integrated optimization method for grouting pressure and flow rate based on diffusion inversion and gas classification as described in claim 1, characterized in that, The process of establishing the force, seepage, and gas coupling control equations and constitutive relations for grouting, surrounding rock, and gas based on the initial parameter field, and setting initial and boundary conditions to obtain an initial multi-field coupled prior model includes: The initial parameters are spatially mapped and interpolated on a unified grid to complete and smooth equivalent permeability, porosity, surrounding rock mechanical properties, fluid viscosity, rheological intensity and interface leakage capacity, thus obtaining a consistent parameter field. The constitutive relationships of grouting, surrounding rock, and gas are determined based on the aforementioned consistency parameter field; A set of governing equations is generated based on the constitutive relation; Initial and boundary conditions are set and the governing equations are coupled to establish an initial multi-field coupled prior model. The initial conditions include: the initial values ​​of the effective stress coupling coefficient, the initial values ​​of the liquid phase pressure, the initial values ​​of the liquid phase saturation, the initial values ​​of the equivalent permeability, the initial values ​​of the liquid phase equivalent viscosity, the initial values ​​of the gas phase saturation, and the initial values ​​of the gas phase pressure. Boundary conditions include: grouting orifice boundary excitation, surrounding rock mechanical boundary, seepage and gas boundary, interface and fracture boundary, compatibility and safety constraints, and monitoring and feedback interface. The grouting orifice boundary excitation is parameterized by boundary excitation parameters. The governing equations include: Mechanical control equations, seepage control equations, and gas control equations; The expression for the governing equations is as follows: ； The expression for the seepage control equation is: ； The expression for the gas control equation is: ； The expression for the initial multi-field coupled prior model is: ； in, For the initial multi-field coupled prior model, It is a displacement vector. The initial field for the displacement vector; The gradient is symmetric for the displacement; To follow the curing indication status The anisotropic stiffness tensor of evolution, The initial values ​​for the anisotropic stiffness tensor; The effective stress coupling coefficient affected by curing. The initial value for the effective stress coupling coefficient; and These are the initial values ​​of the liquid phase pressure and the liquid phase pressure, respectively. It is a unit second-order tensor; Body density; Porosity as it evolves with curing, and This represents the initial value for porosity; and This represents the liquid phase saturation and its initial value. The density of the liquid phase; This refers to the Darcy velocity in the liquid phase. For liquid phase source sinks; The equivalent permeability is modulated by both curing and anisotropy. This is the initial value for the equivalent permeability; and These are the relative permeability functions of the liquid phase and the gas phase, respectively. The equivalent viscosity of the liquid phase as it evolves with solidification. This is the initial value of the liquid phase equivalent viscosity; The liquid phase pressure gradient is g; g is the acceleration due to gravity. and This represents the gas phase saturation and its initial value. This refers to the gas phase density. The equivalent density of the solid framework or adsorption matrix; It is a monotonic function of the adsorption-desorption equilibrium loading that depends on the gas phase pressure; This refers to the Darcy velocity in the gas phase. For gas phase source and sink terms; This refers to the viscosity of the gas phase. and These are the initial values ​​for the gas phase pressure and the gas phase pressure, respectively. For computational domain, The computational domain boundary is defined by the excitation of the grouting orifice boundary, the surrounding rock mechanics boundary, the seepage and gas boundary, the interface and fracture boundary, the compatibility and safety constraints, and the monitoring and feedback interface. This is an initial constraint set that uniformly configures the excitation of mechanics, the seepage and gas boundaries, and the grouting orifice boundaries. The effective stress coupling coefficient is determined by the formation structure, mechanical and seepage parameters, fracture network, and historical construction records, while the liquid phase saturation is determined by the initial gas pressure and content.

7. The integrated optimization method for grouting pressure and flow rate based on diffusion inversion and gas classification as described in claim 6, characterized in that, A grouting boundary excitation scheme is constructed, and the initial multi-field coupled prior model is parameterized with boundary excitation parameters to obtain the final multi-field coupled prior model, including: On the boundary of the initial multi-field coupled a priori model, the boundaries of grouting orifice, surrounding rock mechanics, seepage outer edge, gas outer edge, and fracture and interface are defined, and their geometric range and normal direction are marked to form a parameterizable set of boundaries. Based on the parameterizable boundary set, programmable basic elements integrating pressure and flow are defined for the grouting orifice boundary, standard basic elements of constant pressure, constant flow, and no flux are defined for the outer edge boundary, and displacement and stress basic elements are defined for the mechanical boundary. The parameterized interface of time history, amplitude and switching strategy is uniformly specified to obtain the boundary element library. Based on the boundary element library, the time history of the orifice pressure and flow rate is constructed and applied to obtain the grouting orifice time history excitation pair; Based on the time history excitation pair of the grouting orifice coupled with the fracture and matrix as a consistent expression of normal flux and pressure difference, the interface flux constraint set is obtained; Based on the interface flux constraint set, compatibility checks and scheduling are performed on the mechanical, seepage and gas boundaries to generate an executable boundary time history plan. The relationship between grouting orifice excitation and interface flux is embedded into the governing equations according to the boundary type, and the consistency of the initial conditions and parameter fields is checked to obtain a closed solution body containing time history boundaries. Based on the closed solution body containing the time history boundary, the operator form transformation of the boundary excitation parameterization is performed to obtain the final multi-field coupled prior model; Among them, the time history excitation pair of the grouting orifice includes: the time history of the grouting orifice pressure and the time history of the grouting orifice flow. The interface flux constraint set includes: the normal leakage flux density of the surrounding rock interface and the interface conductivity coefficient; The expressions for the grouting orifice time-history excitation pair, the interface flux constraint set, and the final multi-field coupled prior model are as follows: ； ； ； in, The time history of the grouting orifice pressure; The time history of the grouting orifice flow rate; and The baseline pressure and baseline flow rate for grouting; The pressure and increase in traffic during the phased leap; The pressure and flow rate amplitudes of the pulsating segment; The pulsation frequency; This represents the phase difference between pressure and flow rate. This represents the start time of the corresponding stage; The duration of the phase; A unit rectangular window function is used to define the phase window; denoted as the normal leakage flux density at the surrounding rock interface; The interface conductivity coefficient; The equivalent pore pressure in the interface neighborhood; This is the uncured flux base value; The directional conductivity coefficient in the fully cured state; The unit vector is the interface normal. It is a boundary operator with time and parameter vector as independent variables; For the boundary excitation parameter vector, This is the final multi-field coupled prior model after parameterization of boundary excitations; This is the computational solution body after embedding time history boundaries and completing consistency verification; A set for boundary classification and role labeling; For boundary primitive library; This is the coordinated time history set of grouting orifice pressure and flow rate; For the set of interface flux constraints; For multi-boundary compatible scheduling plans.

8. The integrated optimization method for grouting pressure and flow rate based on diffusion inversion and gas classification as described in claim 7, characterized in that, Based on the aforementioned multi-field coupled prior model, online acquisition of field data and data preprocessing are performed to obtain a time-series dataset, including: Under a unified grid and geometric framework, the state variables and boundary positions of the pressure, flow rate, displacement, gas concentration and pore pressure monitoring channels are mapped one by one with the multi-field coupled prior model to obtain a structured online observation layout; Based on a structured online observation layout, raw data streams are continuously collected from each monitoring channel according to a predetermined sampling period and aligned and deduplicated using a unified time reference to obtain time-synchronized raw time-series records. The original time series records are subjected to quality inspection and saturation, jump, disconnection and drift anomalies are removed to obtain clean time series data; The cleaning time series data is subjected to zero-phase finite bandwidth filtering according to the physical bandwidth of the monitoring indicators and power frequency and mechanical harmonic interference is suppressed to obtain pre-filtered time series data; The pre-filtered time series data of each channel are resampled at a uniform target sampling rate, and short-time missing measurements are filled by unbiased interpolation based on local trends to obtain aligned time series data. The preprocessed synchronization segments from each monitoring channel are merged into a structured time-series dataset in chronological order.

9. The integrated optimization method for grouting pressure and flow rate based on diffusion inversion and gas classification as described in claim 8, characterized in that, Using the final multi-field coupled prior model as the forward modeling skeleton and the time-series dataset as constraints, a joint objective function is constructed. Regularization and physical priors are introduced into this joint objective function to jointly estimate spatiotemporally variable equivalent permeability, porosity, rheological parameters, interfacial leakage coefficient, and diffusion front geometry. Based on the joint estimation results, the posterior distribution of the parameters is determined, including: Under the framework of multi-field coupled forward modeling, an extended state vector containing field states and parameters to be estimated is constructed and time discretized and grid aligned to support the recursive solution of filtering, so as to obtain the extended state definition for assimilation. The parameters to be estimated include spatiotemporally variable equivalent permeability, porosity, rheological parameters, interface leakage coefficient and diffusion front geometry. Based on the time series dataset and monitoring layout, observation operators are defined to map the model state to the observation space, resulting in an observation vector flow and observation operators that match the extended state. The observation vector flow, observation operator, and spatial smoothing regularization term are used to construct a joint objective function and embed physical priors to estimate the parameters to be estimated, thereby obtaining the posterior estimates of the state and parameters at the current time. Based on the posterior estimates of the current state and parameters, the forward modeling skeleton is used to advance the extended state and covariance, and tangentization is used to process nonlinearities to obtain the prior distribution for the next time step. The prior distribution is corrected by observation assimilation and the state and parameters are jointly updated to obtain the posterior distribution at the next time step, and the updated posterior mean and covariance of the state and parameters are obtained. After each update, perform regularized consistent projection and physical domain correction to satisfy nonnegativity, upper bound, anisotropic symmetry and geometric consistency with the diffusion front, and obtain the posterior distribution of parameters; The extended state definition expression, the joint objective function expression, and the expression for the posterior parameter distribution are respectively: ； ； ； in, A timeline for boundaries and external stimuli; This is a preprocessed structured time-series dataset; For the first The extended state vector contains the field state. With the parameters to be estimated It is a nonlinear state transition operator; External stimulus and boundary input; The process noise has a covariance of For observation vectors; For observation operators; To observe the noise, the covariance is... The joint objective function; This is the parameter vector from the previous time step; For spatial smooth covariance; The physical residual vector; For physical residual covariance; and These are the predicted and updated state estimates, respectively; and These are the predicted and updated covariances, respectively. The tangent matrix or set statistics are used to replace the state transition; Tangentization matrix or set statistics substitution for observation operators; Kalman gain; It is an identity matrix.

Images