Four-dimensional safety and environmental protection boundary prediction method, device and equipment for resource development

Abstract

The application provides a four-dimensional safe and environmentally-friendly boundary prediction method, device and equipment for resource development, and relates to the technical field of digital twinning. The method comprises the following steps: obtaining materialization characteristic parameters corresponding to each resource type in a target development scheme, and constructing the materialization characteristic parameters into a unified materialization vector; determining a process parameter tensor according to a mining process operation parameter sequence of the target development scheme; taking the unified materialization vector and the process parameter tensor as driving forces, and coupling calculation of time-space evolution results of equivalent plastic strain, pollutant concentration and vegetation stress index through a pre-established rock mass deformation-fracture sub-model, a groundwater flow-pollutant transport sub-model and an ecological carrier degradation sub-model; and generating a four-dimensional safe and environmentally-friendly time-space boundary based on the time-space evolution results. The application dynamically couples geological disasters, groundwater environmental protection and ecological degradation in terms of safe and environmentally-friendly constraints, and provides a more accurate and quantifiable safe and environmentally-friendly boundary for resource development.

Description

Technical Field

[0001] This invention relates to the field of digital twin technology, and in particular to a method, apparatus and equipment for predicting the four-dimensional safety and environmental protection boundaries of resource development. Background Technology

[0002] Resource development (including coal mining, metal mining, oil and gas extraction, and geothermal development) faces multiple safety and environmental challenges, such as geological disasters, groundwater pollution, and surface ecological degradation. Especially in the development of deep mineral resources, the geological environment is complex, and the mining process often takes place in extreme environments with multiple coupled fields, including high stress, high temperature, high osmotic pressure, and complex chemical fields. With increasing mining depth and intensity, the risks of safety and environmental accidents, such as rock mass instability, water inrush, gas outbursts, surface subsidence, and groundwater pollution, caused by stress redistribution, evolution of rock mass pore pressure, and changes in temperature fields, increase significantly. How to scientifically and dynamically delineate the safe permissible range and environmentally acceptable range for deep resource development activities, and determine the four-dimensional safety and environmental protection boundaries, has become a key technical challenge that urgently needs to be addressed in the field of deep resource safety and green development.

[0003] Currently, most mine safety and environmental protection technical solutions rely on single indicators such as groundwater level monitoring values, microseismic event rates, or surface deformation rates exceeding static thresholds as alarm basis. Essentially, these are post-event / near-event warnings, failing to proactively define the four-dimensional spatiotemporal boundary surfaces—when, at what depth, and within what extent mining is permissible—during the resource development planning and design phase. Furthermore, these methods cannot reflect the temporal evolution characteristics of rock mass mechanical properties and fluid transport patterns under deep multi-field coupling conditions, making it difficult to establish comprehensive criteria for multi-indicator synergy. Consequently, the environmental compliance evaluation of development plans lacks quantifiable and verifiable spatial evidence. Summary of the Invention

[0004] This invention provides a method, apparatus, and equipment for predicting the four-dimensional safety and environmental protection boundaries of resource development, in order to solve the problem that existing technologies cannot accurately predict the safety and environmental protection boundaries of resource development in the case of independent calculation of physical fields and interaction of multiple physical fields.

[0005] In a first aspect, embodiments of the present invention provide a four-dimensional safety and environmental protection boundary prediction method for resource development, including: Obtain the materialized characteristic parameters corresponding to each resource type in the target development plan and construct them into a unified materialized vector; Based on the sequence of operational parameters for the mining process of the target development scheme, determine the process parameter tensor; Using the unified materialization vector and the process parameter tensor as driving forces, the spatiotemporal evolution results of equivalent plastic strain, pollutant concentration and vegetation stress index are coupled and calculated through the pre-established rock mass deformation-fracture sub-model, groundwater flow-pollutant transport sub-model and ecological carrier degradation sub-model. Based on the spatiotemporal evolution results, a four-dimensional safe and environmentally friendly spatiotemporal boundary is generated.

[0006] In one possible implementation, the unified materialization vector and the process parameter tensor are used as driving forces to couple and calculate the spatiotemporal evolution results of equivalent plastic strain, pollutant concentration, and vegetation stress index through a pre-established rock mass deformation-fracture sub-model, a groundwater flow-pollutant transport sub-model, and an ecological carrier degradation sub-model, including: For each spatial discrete grid, the unified materialization vector and the process parameter tensor are input into the rock mass deformation-fracture sub-model, and the water-conducting fracture zone height, equivalent plastic strain and fracture permeability of the spatial discrete grid are output. The fracture permeability is updated to the groundwater flow-contaminant transport sub-model, and the unified materialization vector and the process parameter tensor are input into the updated groundwater flow-contaminant transport sub-model. The groundwater drawdown and rock pore pressure of the spatial discrete grid are output. When the height of the water-conducting fracture zone is greater than the fracture threshold, the contaminant concentration of the spatial discrete grid is output. The groundwater drawdown is updated to the ecological carrier degradation sub-model. The unified materialization vector and the process parameter tensor are input into the updated ecological carrier degradation sub-model, and the vegetation stress index of the spatial discrete grid is output. The rock mass pore pressure is updated to the rock mass deformation-fracture sub-model, and the iterative calculation for the next time point is initiated. Based on the time series values ​​of the equivalent plastic strain, pollutant concentration, and vegetation stress index corresponding to each spatial discrete grid, the spatiotemporal evolution results are determined.

[0007] In one possible implementation, generating a four-dimensional safe and environmentally friendly spatiotemporal boundary based on the spatiotemporal evolution result includes: Based on the relationship between the equivalent plastic strain and the critical subsidence value within each spatial discrete grid, a time-varying geological hazard mining limit Boolean field is generated. Based on the relationship between the pollutant concentration in each spatial discrete grid and the preset pollution threshold, a groundwater environmental Boolean field that changes over time is generated. Based on the relationship between the vegetation stress index and the preset vegetation stress threshold in each spatial discrete grid, a three-dimensional ecological degradation dynamic Boolean field that changes over time is generated. Based on the geological hazard restricted mining Boolean field, the groundwater environmental protection Boolean field, and the three-dimensional ecological degradation restricted dynamic Boolean field, the four-dimensional safe and environmentally friendly spatiotemporal boundary is generated.

[0008] In one possible implementation, generating the four-dimensional safe and environmentally friendly spatiotemporal boundary based on the geological hazard-limited mining Boolean field, the groundwater environmental protection Boolean field, and the three-dimensional ecological degradation-limited dynamic Boolean field includes: By performing a logical OR operation on the geological hazard restricted mining Boolean field, the groundwater environmental protection Boolean field, and the three-dimensional ecological degradation restricted movement Boolean field, a comprehensive prohibited entry Boolean field is obtained; The four-dimensional safe and environmentally friendly spatiotemporal boundary is determined based on the complement of the comprehensive forbidden Boolean field.

[0009] In one possible implementation, the formula for the rock mass deformation-fracture sub-model includes: ∇ ·s + p g b = 0 ; s = ( 1 - D )· C : e e ; Ḋ = f ( e eq , s 1 , s 3 , M, H f ); Hf =∑ n i=1 hi ⋅ I ( Yes < Scr ); k f = k 0 · exp ( a · e v ) ( e v > e cr ); in, s For rock mass stress tensor, r g The rock mass density is driven by the unified materialization vector. b For the volume forces of the rock mass, ∇ · For divergence operators,D For rock mass damage variables, C Let be the rock mass elastic stiffness tensor. e e Let be the elastic strain tensor of the rock mass, and be the double dot product. e eq For the equivalent plastic strain of the rock mass, s 1 The maximum principal stress of the rock mass, s 3 The minimum principal stress of the rock mass, M The sampling height is driven by the process parameter tensor. H f The height of the water-conducting fracture zone, hi For the first i Thickness of rock layers, Yes For the first i The value of the drop in rock strata, Scr This represents the critical subsidence value of the rock mass. I For indicator functions, n This represents the total number of overburden layers. k f The fracture permeability, k 0 The initial permeability of the rock mass. α The damage-permeability coupling coefficient is... e v For the volumetric strain of the rock mass, exp ( · ) is an exponential function, e cr The critical volumetric strain of the rock mass is given.

[0010] In one possible implementation, the formula for the groundwater flow-contaminant transport sub-model includes: ∂ ( θp ) / ∂t = ∇⋅( k f / m⋅ ∇ p ) +Q s ; v = -k f / m; (∇ p + p y g ∇ z ); ∂ ( θC ) / ∂t = ∇·( θD· ∇ C ) - ∇( θvC ) + θR ( C ) + q s C s ; in, i The rock mass porosity is driven by the unified materialization vector. p The pore pressure of the rock mass is... k f The fracture permeability, m For the unified materialization vector-driven hydrodynamic viscosity, Q s Let ∇ be the source and sink terms of the seepage equation driven by the unified materialization vector, ∇ be the gradient operator, and ∇⋅ be the divergence operator. v The average flow velocity of the fluid in the void. r y For the fluid density driven by the unified materialization vector, g It is the acceleration due to gravity. z For elevation, C The concentration of the pollutant, D The hydrodynamic dispersion coefficient tensor driven by the unified materialization vector. R ( C () represents a chemical reaction term. C s The source and sink concentrations are driven by the process parameter tensor. q s These are the source and sink terms of the concentration equation driven by the unified materialization vector.

[0011] In one possible implementation, the formula for the ecological carrier degradation sub-model includes: E stress ( x,y,t ) =w 1 ·f 1 (Δ h / h max ) + w 2 ·f 2 (Δ i r / i fc ) +w 3 ·f 3 ( PET / P ); in, Estress ( x,y,t ) for time t Plane coordinates ( x , y The vegetation stress index described below, Dh For the aforementioned groundwater drawdown, h max The maximum root water absorption depth of vegetation driven by the unified materialization vector. i r For the residual content, i fc For the unified materialized vector-driven field water holding capacity, PET For potential evaporation, P For precipitation, f 1 The first vegetation type response coefficient is driven by the unified materialization vector. f 2 The second vegetation type response coefficient is driven by the unified materialization vector. f 3 The response coefficient of the third vegetation type driven by the unified materialization vector. w 1 The first weight coefficient driven by the unified materialization vector. w 2 The second weighting coefficient is driven by the unified materialization vector. w 3 The third weight coefficient driven by the unified materialization vector satisfies w 1 + w 2 + w 3 =1.

[0012] In one possible implementation, obtaining the materialized characteristic parameters corresponding to each resource type in the target development scheme and constructing them into a unified materialized vector includes: Obtain the physicochemical property parameters corresponding to the solid phase resources, liquid phase resources and gas phase resources in the target development scheme; The physicochemical property parameters corresponding to the solid phase resources, liquid phase resources and gas phase resources are mapped to a unified physical property coordinate system for standardization to obtain standardized solid phase parameters, standardized liquid phase parameters and standardized gas phase parameters. The standardized solid-phase parameters, the standardized liquid-phase parameters, and the standardized gas-phase parameters are concatenated and multiplied by the transformation matrix to obtain the unified physicalization vector.

[0013] Secondly, embodiments of the present invention provide a four-dimensional safety and environmental protection boundary prediction device for resource development, comprising: The unified dimension conversion module is used to obtain the materialized characteristic parameters corresponding to each resource type in the target development plan and construct them into a unified materialized vector. The process parameter acquisition module is used to determine the process parameter tensor based on the sequence of mining process operation parameters of the target development scheme; The predicted spatiotemporal evolution acquisition module is used to use the unified materialization vector and the process parameter tensor as driving forces, and through the pre-established rock mass deformation-fracture sub-model, groundwater flow-pollutant transport sub-model and ecological carrier degradation sub-model, to couple and calculate the spatiotemporal evolution results of equivalent plastic strain, pollutant concentration and vegetation stress index. The safe and environmentally friendly spatiotemporal boundary generation module is used to generate a four-dimensional safe and environmentally friendly spatiotemporal boundary based on the spatiotemporal evolution results.

[0014] Thirdly, embodiments of the present invention provide a four-dimensional safety and environmental protection boundary prediction device for resource development, including a memory and a processor. The memory stores a computer program, and the processor executes the computer program to implement the four-dimensional safety and environmental protection boundary prediction method for resource development as described in the first aspect or any possible implementation of the first aspect.

[0015] In this embodiment of the invention, by constructing the physical property parameters corresponding to each resource type into a unified physical vector, various parameters of solid, liquid, and gas with vastly different dimensions and magnitudes are mapped to a unified physical property coordinate system. This allows subsequent models to obtain corresponding data from the same data interface, supporting generalized prediction of multiphase resources under the same algorithm framework. By using the unified physical vector and process parameter tensor as driving forces, the spatiotemporal evolution results of equivalent plastic strain, pollutant concentration, and vegetation stress index are coupled and calculated on the rock mass deformation-fracture sub-model, the groundwater flow-pollutant transport sub-model, and the ecological carrier degradation sub-model. This solves the technical problem in the prior art where multiple physical fields are calculated independently and cannot dynamically reflect the interactive influence of multiple physical fields, providing a precise and quantifiable four-dimensional spatiotemporal prediction basis for resource development. Attached Figure Description

[0016] Figure 1 This is a flowchart illustrating the implementation of the four-dimensional safety and environmental protection boundary prediction method for resource development provided in this embodiment of the invention. Figure 2 This is a schematic diagram of the structure of the four-dimensional safety and environmental protection boundary prediction device for resource development provided in this embodiment of the invention; Figure 3 This is a schematic diagram of a four-dimensional safety and environmental protection boundary prediction device for resource development provided in an embodiment of the present invention. Detailed Implementation

[0017] The embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0018] like Figure 1 As shown, the four-dimensional safety and environmental protection boundary prediction method for resource development provided in this embodiment of the invention includes: Step 101: Obtain the materialized characteristic parameters corresponding to each resource type in the target development scheme and construct them into a unified materialized vector.

[0019] In one embodiment, step 101 specifically includes: Obtain the physicochemical properties of solid, liquid, and gaseous resources in the target development plan.

[0020] The physicochemical property parameters corresponding to solid, liquid, and gaseous resources are mapped to a unified physical property coordinate system for standardization, resulting in standardized solid parameters, standardized liquid parameters, and standardized gaseous parameters.

[0021] The standardized solid-phase parameters, standardized liquid-phase parameters, and standardized gas-phase parameters are concatenated and multiplied by the transformation matrix to obtain a unified physicochemical vector.

[0022] This step standardizes the physical and chemical parameters of solid-liquid-gas heterogeneous systems, enabling subsequent models to obtain relevant data from the same data interface and supporting generalized prediction of multiphase resources under the same algorithm framework.

[0023] In this embodiment, the physicochemical property parameters corresponding to solid-phase resources, liquid-phase resources, and gas-phase resources are mapped to a unified physical property coordinate system for standardization, resulting in standardized solid-phase parameters, standardized liquid-phase parameters, and standardized gas-phase parameters. Specifically, this includes: Solid phase resource physicochemical property parameters are constructed into solid phase resource components. The solid phase resource components are then divided element-wise by the reference rock mechanics parameter vector to obtain standardized solid phase parameters.

[0024] The physicochemical properties of liquid phase resources are constructed as liquid phase resource components. The difference between the liquid phase resource components and the standard state condition vector is calculated. The difference is then divided by the standard state condition vector element by element to obtain the standardized liquid phase parameters.

[0025] The physical and chemical properties of gaseous resources are constructed into gaseous resource components. The gaseous resource components are then divided element-wise by the critical state parameter vector to obtain standardized gaseous parameters.

[0026] For example, in step 101, the resource type identifier of the target development scheme is received, the physicochemical characteristic parameters of the corresponding solid, liquid, and gaseous resources are retrieved, and a unified physicochemical characteristic vector is constructed. V physchem And map it to a unified physical property coordinate system.

[0027] Among them, the physicochemical properties of solid phase resources include parameters of coal / metal ore: density, hardness, tensile / compressive strength of roof and floor, sulfide content, and acid potential (AP).

[0028] Physicochemical properties of liquid resources include parameters related to oil / gas / groundwater: viscosity, saturation pressure, CH4 solubility, Cl⁻ / SO4²⁻ / heavy metal concentration, and corrosive CO2 partial pressure.

[0029] The physicochemical properties of gaseous resources include relevant parameters of natural gas / CH4: compressibility factor, diffusion coefficient, adsorption desorption isotherm, and prominence tendency index.

[0030] Correspondingly, the solid phase resource components (coal / metal ore) are: V solid =[ r , H , s t , s c , S py , AP ]; in, r The density of solid-phase resources (kg / m³) H The Mohs hardness of solid-phase resources. s t The tensile strength (MPa) of the top and bottom plates of the solid phase resource. s c The compressive strength (MPa) of solid-phase resources, S py The sulfide content (wt%) of solid-phase resources. AP The potential for acid production from solid-phase resources (kg H2SO4 / t).

[0031] The liquid phase resource components (oil and gas / groundwater) are: V liquid = [ m , P sat , S CH4 , C Cl , C SO4 , C HM , P CO2 ]; in, m The dynamic viscosity (Pa·s) of the liquid phase resource.P sat The saturation pressure (MPa) of the liquid phase resource. S CH4 CH4 solubility (mg / L) C Cl The concentration of Cl⁻ (mg / L) C SO4 SO4²⁻ concentration (mg / L) C HM The concentration of heavy metals (mg / L) P CO2 This represents the corrosive CO2 partial pressure (kPa).

[0032] The gas phase resource components (natural gas / coalbed methane) are: V gas = [ Z, D, L a , I out ]; in, Z The compressibility factor for gaseous resources. D Let be the diffusion coefficient (m² / s) of the gas phase resource. L a The parameters of the adsorption analysis isotherm for gas phase resources (Langmuir parameters). I out It is an index indicating a prominent tendency for gas phase resources.

[0033] Accordingly, the standardized solid-state parameters are: V solid ’ = V solid / V ref_solid ; in, V ref_solid This is a reference rock mechanics parameter vector.

[0034] Standardized liquid chromatography parameters are: V liquid ’ = ( V liquid - V std ) / V std ; in, V std This is the standard state condition vector.

[0035] The standardized gas phase parameters are: V gas ’ = V gas / V critical ; in, V critical This is the critical state parameter vector.

[0036] The unified materialization vector is: V physchem = M eq · [ V solid ’ ; V liquid ’ ; V gas ’ ]; in, M eq The transformation matrix is ​​for normalization.

[0037] Step 102: Determine the process parameter tensor based on the sequence of operation parameters of the mining process in the target development plan.

[0038] In one embodiment, step 102 specifically includes: Based on the target development scheme, the sequence of operating parameters for the mining process is analyzed, including: the sequence of operating parameters for solid phase mining, liquid phase mining, gas phase mining, and geothermal mining.

[0039] The process parameter tensors are assembled into a time dimension based on the operating parameter sequences of solid phase mining, liquid phase mining, gas phase mining, and geothermal mining processes.

[0040] For example, in step 102, the sequence of process operation parameters is parsed from the development engineering design document and assembled into a process parameter tensor. T process ( t ).

[0041] The operating parameters for solid phase mining include: mining height, advance speed, amount of residual coal in the goaf, and filling rate.

[0042] The operating parameters for the liquid phase extraction process include: injection pressure, flow rate, fracturing fluid chemical composition, and reinjection temperature.

[0043] Operating parameters for gas phase extraction include: desorption rate, extraction negative pressure, and wellhead back pressure.

[0044] Operating parameters for geothermal extraction include: circulation flow rate, temperature difference between injection and production wells, and operating life.

[0045] Accordingly, the sequence of operating parameters for the solid phase mining process is as follows: T mining ( t )= [ M ( t ), v adv ( t ), R goaf ( t ), or fill ( t )]; in, M ( t () represents the mining height (m). v adv ( t () represents the propulsion speed (m / d). R goaf ( t The figure represents the amount of coal remaining in the goaf (t). or fill ( t ) represents the filling rate (%).

[0046] The sequence of operating parameters for the liquid phase extraction process is as follows: T liquid ( t ) = [ P inj ( t ), Q inj ( t ), C frac ( t ), T reinj ( t )]; in: P inj ( t The injection pressure is (MPa). Q inj ( t () represents the injection displacement (m³ / d). C frac ( t () represents the chemical composition vector of the fracturing fluid.T reinj ( t ( ) represents the recharge temperature (°C).

[0047] The sequence of operating parameters for gas phase extraction is as follows: T gas ( t ) = [ v des ( t ), P vac ( t ), P back ( t )]; in, v des ( t The desorption rate is (m³ / d). P vac ( t The extraction negative pressure (kPa) is the pressure at which the pump is applied. P back ( t ) represents the wellhead back pressure (MPa).

[0048] The sequence of operating parameters for geothermal extraction is as follows: T geo ( t ) = [ v circ ( t ),Δ T wp ( t ), Y op ( t )]; in, v circ ( t ) represents the circulation velocity (m³ / h), Δ T wp ( t ( ) represents the temperature difference between the injection and production wells (°C). Y op ( t ) represents the operating years (a).

[0049] Step 103: Using the unified materialization vector and process parameter tensor as the driving force, the spatiotemporal evolution results of equivalent plastic strain, pollutant concentration and vegetation stress index are coupled and calculated through the pre-established rock mass deformation-fracture sub-model, groundwater flow-pollutant transport sub-model and ecological carrier degradation sub-model.

[0050] In one embodiment, step 103 specifically includes: For each spatial discrete grid, the unified materialization vector and process parameter tensor are input into the rock mass deformation-fracture sub-model, and the height of the water-conducting fracture zone, equivalent plastic strain, and fracture permeability of the spatial discrete grid are output.

[0051] The fracture permeability is updated to the groundwater flow-contaminant transport sub-model, and the unified materialization vector and process parameter tensor are input into the updated groundwater flow-contaminant transport sub-model. The groundwater drawdown and rock pore pressure of the spatial discrete grid are output. When the height of the water-conducting fracture zone is greater than the fracture threshold, the contaminant concentration of the spatial discrete grid is output.

[0052] The groundwater drawdown is updated to the ecological carrier degradation sub-model. The unified materialization vector and process parameter tensor are input into the updated ecological carrier degradation sub-model, and the vegetation stress index of the spatial discrete grid is output.

[0053] The rock mass pore pressure is updated to the rock mass deformation-fracture sub-model, and iterative calculations are performed at the next time point. The spatiotemporal evolution results are determined based on the time series values ​​of equivalent plastic strain, pollutant concentration, and vegetation stress index corresponding to each spatial discrete grid.

[0054] This step, through the step-by-step updating of three sub-models, achieves the synchronous evolution prediction of equivalent plastic strain, pollutant transport, and vegetation stress under mining conditions, ensuring data consistency of different physical fields on the same spatiotemporal scale, and providing technical support for fully coupled numerical simulation of groundwater-ecology co-evolution in mining areas.

[0055] In this embodiment, the formula for the rock mass deformation-fracture sub-model includes: ∇ ·s + p g b = 0 ; s = ( 1 - D )· C : e e ; Ḋ = f ( e eq , s 1 , s 3 , M, H f ); Hf =∑ n i=1 hi ⋅ I ( Yes < Scr ); k f = k 0 · exp ( a · e v ) ( e v > e cr ); in, s For rock mass stress tensor, r g To unify the rock mass density driven by the materialization vector, b For the volume forces of the rock mass, ∇ · For divergence operators, D For rock mass damage variables, C Let be the rock mass elastic stiffness tensor. e e Let be the elastic strain tensor of the rock mass, and be the double dot product. e eq For the equivalent plastic strain of the rock mass, s 1 The maximum principal stress of the rock mass, s 3 The minimum principal stress of the rock mass, M For process parameter tensor-driven sampling height, H f The height of the water-conducting fracture zone, hi For the first i Thickness of rock layers, Yes For the first i The value of the drop in rock strata, Scr This represents the critical subsidence value of the rock mass. I For indicator functions, n This represents the total number of overburden layers. k f For fracture permeability, k 0 The initial permeability of the rock mass (m²) α The damage-permeability coupling coefficient is... e v For the volumetric strain of the rock mass, exp ( · ) is an exponential function, e cr The critical volumetric strain of the rock mass is given.

[0056] The rock mass deformation-fracture sub-model can accurately calculate the height of the water-conducting fracture zone, the equivalent plastic strain, and the fracture permeability, providing a precise guarantee for boundary generation and the updating of the groundwater flow-contaminant transport sub-model.

[0057] For example, the height of the water-conducting fracture zone can also be obtained by directly reading fracture network evolution data through a commercial software interface and extracting the maximum development height of vertically penetrating fractures; the equivalent plastic strain is obtained by internal calculation in the rock mass deformation-fracture sub-model through numerical simulation software.

[0058] In this embodiment, the formula for the groundwater flow-contaminant transport sub-model includes: ∂ ( θp ) / ∂t = ∇⋅( k f / m⋅ ∇ p ) +Q s ;; v = -k f / m; (∇ p + p y g ∇ z ); ∂ ( θC ) / ∂t = ∇·( θD· ∇ C ) - ∇( θvC ) + θR ( C ) + q s C s ; in, i To unify the materialized vector-driven rock mass porosity, p This refers to the pore pressure of the rock mass. k f For fracture permeability, m To unify the materialized vector-driven fluid dynamic viscosity, Q s To unify the source and sink terms (1 / s) of the seepage equation driven by materialization vectors, ∇ is the gradient operator, and ∇⋅ is the divergence operator. v The average flow velocity of the fluid in the void. r y To unify the materialization vector-driven fluid density, g It is the acceleration due to gravity. z For elevation, C The concentration of pollutants (mg / L) D To unify the hydrodynamic dispersion coefficient tensor (m² / s) driven by the materialization vector, R ( C() represents the chemical reaction term (mg / L / s). C s Source and sink concentrations (mg / L) driven by process parameter tensors. q s The source and sink terms (1 / s) of the concentration equation are unified by the materialization vector.

[0059] The groundwater flow-contaminant transport sub-model can accurately calculate groundwater drawdown, rock pore pressure, and contaminant concentration, providing precise assurance for updating the boundary generation and rock deformation-fracture sub-model and the ecological carrier degradation sub-model.

[0060] For example, the groundwater flow-contaminant transport sub-model also includes the following formula: i = i m + i f ; D = D m + α L |v| ; The finite volume method (FVM) is used to solve for pollutant concentrations in the groundwater flow-pollutant transport sub-model on a shared spatial discrete grid. ∫ V ( ∂(θC) / ∂t ) dV =∮ Γ θD· ∇ C n 1 dΓ-∮ Γ θCv · n 1 dΓ + ∫ V θR(C) dV + ∫ V q s C s dV ; in, V For groundwater volume, C This is the boundary of the groundwater volume. n 1 The boundary normal vector.

[0061] Discretize the time using implicit Euler or Crank-Nicolson schemes: ( i {t+1} C t+1 - i t C t ) / Δ t = L ( C t+1 ) + S t+1 ; in, t At the current moment, Δ t For time step.

[0062] In this embodiment, the formula for the ecological carrier degradation sub-model includes: E stress ( x,y,t ) =w 1 ·f 1 (Δ h / h max ) + w 2 ·f 2 (Δ i r / i fc ) +w 3 ·f 3 ( PET / P ); in, E stress ( x,y,t ) for time t Plane coordinates ( x , y The vegetation stress index under ) Dh To reduce the depth of groundwater, h max To standardize the maximum root water absorption depth of vegetation driven by materialized vectors, i r For the residual content, i fc To unify the materialized vector-driven field water holding capacity, PET Potential evapotranspiration (mm / d) was calculated using the Penman-Monteith formula. P Rainfall (mm / d) f 1To unify the response coefficients of the first vegetation type driven by the materialization vector, f 2 To unify the response coefficients of the second vegetation type driven by the materialization vector, f 3 To unify the response coefficients of the third vegetation type driven by materialization vectors, f 1 , f 2 , f 3 Based on the mapping relationship (linear, logistic, or threshold function) of vegetation type response function, w 1 The first weight coefficient is driven by the unified materialized vector. w 2 The second weighting coefficient is driven by the unified materialization vector. w 3 To unify the materialized vector-driven third weight coefficient, satisfying w 1 + w 2 + w 3 =1.

[0063] The ecological carrier degradation sub-model can accurately calculate the vegetation stress index, providing precise assurance for boundary generation.

[0064] For example, the ecological carrier degradation sub-model also includes the following formula: Δ h ( x , y , t ) = Q / (4π T ) · W ( u ); u = r ² S / (4 Tt ); Where, Δ h ( x , y , t This refers to the spatiotemporal evolution of groundwater drawdown. Q Pumping volume (m³ / d) T The hydraulic conductivity is (m² / d). W ( u ) is the Theis well function. r Distance from the pumping well (m) S This represents the water storage coefficient.

[0065] For multi-well interference or complex boundary conditions, numerical solutions are used: S s ∂ h / ∂ t = ∇·( T ·∇ h ) + Q s ; in, S s The water storage rate is 1 / m.

[0066] Δ i r ( x , y , t ) = i fc i ( z root , t ); Where, Δ i r ( x , y , t (This refers to the spatiotemporal evolution of the residuals, which includes quantitative changes.) i fc It refers to field water holding capacity (volume water content). c i ( z root , t ) represents the root region depth z root Volumetric moisture content at that location.

[0067] Calculated from soil moisture characteristic curve and groundwater level depth: i ( z root , t ) = i r + ( i s -θ r) / [1 + | a·ψ ( z root , t )| n 2 ] 1-1 / n 2 ; ψ ( z root ,t ) = -Δ h ( x , y , t )+ z root ; in, i s Where α is the saturated water content, n2 is the inverse parameter of the intake suction force, and n2 is the porosity distribution index. ψ ( z root , t ) represents the matrix potential (m).

[0068] Step 104: Based on the spatiotemporal evolution results, generate a four-dimensional safe and environmentally friendly spatiotemporal boundary.

[0069] In one embodiment, step 104 specifically includes: Based on the relationship between the equivalent plastic strain and the critical subsidence value within each spatial discrete grid, a time-varying geological hazard mining limit Boolean field is generated. Based on the relationship between the pollutant concentration in each spatial discrete grid and the preset pollution threshold, a groundwater environmental Boolean field that changes over time is generated. Based on the relationship between the vegetation stress index and the preset vegetation stress threshold in each spatial discrete grid, a three-dimensional ecological degradation dynamic Boolean field that changes over time is generated. Based on the geological disaster-limited mining Boolean field, the groundwater environmental protection Boolean field, and the three-dimensional ecological degradation-limited dynamic Boolean field, a four-dimensional safe and environmentally friendly spatiotemporal boundary is generated.

[0070] This step constructs a dynamic Boolean field encompassing three dimensions: geological hazards, groundwater environmental protection, and ecological degradation. This allows for consideration of safety and environmental constraints in these three areas, providing quantifiable boundaries for resource development.

[0071] In this embodiment, a four-dimensional safe and environmentally friendly spatiotemporal boundary is generated based on the geological hazard-limited mining Boolean field, the groundwater environmental protection Boolean field, and the three-dimensional ecological degradation-limited dynamic Boolean field, specifically including: By performing a logical OR operation on the geological hazard restricted mining Boolean field, the groundwater environmental protection Boolean field, and the three-dimensional ecological degradation restricted movement Boolean field, a comprehensive prohibited entry Boolean field is obtained.

[0072] Based on the complement of the comprehensive prohibited Boolean field, the four-dimensional safe and environmentally friendly spacetime boundary is determined.

[0073] This step achieves integrated quantification and dynamic delineation of multiple temperature-related environmental safety factors through Boolean field fusion and boundary generation.

[0074] For example, in step 104, for each time step t, the critical state Boolean field of the rock mass deformation-fracture sub-model, the groundwater flow-pollutant transport sub-model, and the ecological carrier degradation sub-model is calculated on the shared spatial discrete grid.

[0075] Geological hazards restrict mining in Boer fields: B geo ( x , y , z , t ) = 1, e eq ( x , y , z , t ) ≥ ε cr ; B geo ( x , y , z , t ) = 0, e eq ( x , y , z , t )<ε cr ; in, B geo ( x , y , z , t (This is a geological hazard-restricted mining area.) e eq ( x , y , z , t ) represents the spatiotemporal evolution of equivalent plastic strain, ε cr This is the critical subsidence value.

[0076] Groundwater environmental protection site: B env ( x , y , z , t ) = 1, C ( x , y , z , t ) ≥ C limit ; Benv ( x , y , z , t ) = 0, C ( x , y , z , t )< C limit ; in, B env ( x , y , z , t (This is a groundwater environmental protection site.) C ( x , y , z , t ( ) represents the spatiotemporal evolution of pollutant concentration. C limit To preset the pollution threshold, C limit For water / air pollution, the appropriate standards can be determined based on relevant groundwater quality standards, such as total dissolved solids (TDS) ≤ 1000 mg / L, and heavy metal concentrations according to the limits for each element. For geothermal / gas storage facilities, there is no uniform environmental upper limit for geothermal water temperature, but high-temperature extraction (e.g., >90℃) poses environmental risks. The pressure threshold for gas storage facilities can be determined on a case-by-case basis according to relevant underground gas storage design specifications or compressed air energy storage power station underground gas storage design specifications.

[0077] Ecological degradation and limited activity: Boolean field B eco ( x , y , t ) = 1, E stress ( x , y , t ) ≥ E threshold ; B eco ( x , y , t ) = 0, E stress ( x , y , t )< E threshold ; in,B eco ( x , y , t (This is a limited area for ecological degradation.) E stress ( x , y , t () represents the spatiotemporal evolution of the vegetation stress index. E threshold The vegetation stress threshold. E threshold Determined based on vegetation type and ecological function zoning, such as shrubs in arid areas. E threshold = 0.6, grassland area E threshold = 0.5.

[0078] Projecting a two-dimensional ecological degradation-limited Boolean field into three-dimensional space: B eco3D ( x,y,z,t ) = B eco ( x , y , t ) · I(z ≤ z surf ) ; in, B eco3D ( x,y,z,t (This is a three-dimensional ecological degradation limit Boolean field.) z surf Let I be the ground elevation (m), and I(·) be the indicator function.

[0079] Calculate the comprehensive forbidden Boolean field: B forbid ( x,y,z,t ) = B geo ( x,y,z,t ) ∨B env ( x,y,z,t ) ∨B eco3D ( x,y,z,t ); in, B forbid ( x,y,z,t (This is a comprehensive ban on entry into the Boer field.)

[0080] The safety boundary is the complement of the forbidden Boolean field: Ω safe ( t ) = {( x,y,z ) | B forbid ( x,y,z,t ) = 0}; Among them, Ω safe ( t () is the safety boundary.

[0081] For Ω safe ( t Perform morphological closing operations to eliminate numerical noise: Ω safe ’ ( t = close(Ω) safe ( t ), SE); Among them, Ω safe ’ ( t ) represents the safe boundary after denoising, close(·) is the morphological closing operation (dilation followed by erosion), and SE is the structuring element (such as a 3×3×3 cube).

[0082] The moving cubes algorithm is used to extract isosurfaces and generate boundary triangulation meshes. ∂Ω safe ( t ) = MarchingCubes(Ω safe ’ ( t ), iso value =0.5); Among them, ∂Ω safe ( t () represents the boundary of the safe area, MarchingCubes(·) is the algorithm for moving cubes, and iso value The isosurface threshold.

[0083] Further categorize the restricted areas: Prohibited objects: B geo Areas with a value of 1 are marked as "Mining Absolutely Prohibited".

[0084] Pollution protection body: B env = 1 and B geo Areas with a value of 0 are marked as "Restricted Mining / Requires Protective Measures".

[0085] Ecologically confined organisms: B eco3D = 1 and B geo = Benv Areas with a value of 0 are marked as "Surface Activity Restricted Zones".

[0086] In some embodiments, the four-dimensional safety and environmental protection boundary prediction method for resource development further includes: Step 105: Output and compliant delivery of four-dimensional safety and environmental protection spatiotemporal boundaries.

[0087] Ω safe ( t The data is rasterized into a 3D geological GIS layer, generating a safety and environmental protection boundary report, which includes XYZ coordinates, effective time windows, sensitivity ranges, and a list of key assumptions.

[0088] In one embodiment, step 105 specifically includes: Ω safe ( t According to time step Δ t (e.g., 1 month, 1 quarter) Discretized into time-series raster data. Each raster cell stores a security status code: 0: prohibited from acquisition / entry (red), 1: restricted from acquisition / entry (yellow), 2: permitted from acquisition / entry (green).

[0089] Generate a 3D geological layer in GeoTIFF format, containing the following bands: Band 1: Safety status code (0 / 1 / 2); Band 2: Geological hazard risk level (0-5); Band 3: Groundwater pollution risk level (0-5); Band 4: Ecological degradation risk level (0-5); Band 5: Confidence level (0-100%).

[0090] Reports are automatically generated based on the approval templates from mining / geothermal / oil and gas authorities, including: XYZ coordinate table, effective time window, sensitivity range, list of key assumptions, and standard interface.

[0091] The XYZ coordinate table represents the coordinates of key nodes on the boundary surface (WGS84 or CGCS2000 coordinate system); the effective time window represents the effective start and end time of the safety status of each area; the sensitivity interval represents the displacement range of the boundary surface when the key parameters change by ±10%; the list of key assumptions includes: model assumptions, parameter sources, and uncertainty descriptions; standard interface: supports exporting to GeoTIFF, GeoJSON, KML, and CityGML formats, and is compatible with ArcGIS, QGIS, and SuperMap platforms.

[0092] The following specific embodiments illustrate the four-dimensional safety and environmental protection boundary prediction method for resource development provided by the present invention: Example 1: Prediction of safety and environmental protection boundaries in deep coal mining.

[0093] A certain mining area is planned to mine the No. 3 coal seam, with a burial depth of 800-1200m, a mining height of 3.5m, a working face length of 200m, and an advance speed of 5m / d. The overlying strata contain Quaternary loose aquifers and bedrock fissure aquifers, and the surface is a grassland ecosystem.

[0094] S1, Resource Form - Material Characteristics Encoding.

[0095] Resource type identification: Deep coal mining (mainly solid phase). Physicochemical property parameters retrieved: No. 3 coal seam: ρ = 1350 kg / m³, H = 2.5, σ t = 1.2 MPa, σ c = 18 MPa, S py = 2.1 wt%, AP = 8.5 kg H2SO4 / t; Top sandstone: ρ = 2600 kg / m³, H = 6.5, σ t = 5.8 MPa, σ c = 85 MPa; Base mudstone: ρ = 2400 kg / m³, H = 3.0, σ t = 2.1 MPa, σ c = 35 MPa.

[0096] Build V physchem = [1350, 2.5, 1.2, 18, 2.1, 8.5, 2600, 6.5, 5.8, 85, 2400, 3.0, 2.1, 35].

[0097] S2, Process parameter extraction.

[0098] Extracted from the target development plan: T process ( t ) = [ M ( t =3.5m, v adv ( t =5m / d, R goaf ( t =150t, or fill ( t [0%] (This example is a fully mechanized mining face without backfilling) S3, Multi-field Coupled Simulation.

[0099] S3a rock mass deformation-fracture sub-model: A three-dimensional numerical model was built using FLAC3D. The model size was 1000m×500m×600m (length×width×height). The mesh was divided as follows: 5m×5m×2m near the working surface and 20m×20m×5m in the far field.

[0100] Constitutive model: Mohr-Coulomb elastoplastic model + strain softening damage.

[0101] Boundary conditions: bottom fixed, all four sides constrained by normals, top free.

[0102] Initial stress calculation: s v = ρgz , s h = K 0 · s v , K 0 = 0.8.

[0103] in, s v For vertical stress, z For burial depth, s h For horizontal stress, K 0 This is the lateral pressure coefficient.

[0104] Simulation results: When the working face advances to 500m, the height of the water-conducting fracture zone... Hf = 68m, fracturing ratio 19.4. The plastic zone extends to the surface ( B geo = 1).

[0105] S3b Groundwater Flow-Contaminant Transport Sub-Model: when Hf = When the aquifer bottom plate is buried at a depth of 45m, the pollutant transport simulation is triggered.

[0106] Aquifer parameters: k = 5.2 × 10⁻ 4 m / s, i = 0.25, α L = 10m, D m = 1×10⁻ 9 m² / s.

[0107] in, k Permeability coefficient, α L For longitudinal dispersion, Dm is the molecular diffusion coefficient.

[0108] Pollution source: accumulated water in the mined-out area, TDS = 3500 mg / L, SO4²⁻ = 1200 mg / L, pH = 3.2.

[0109] The FVM method was used for the solution, with a time step of Δt = 1d and a total simulation time of 5 years.

[0110] Results: The leading edge of the pollution plume reached the downstream boundary at 3 years, with a concentration of... Cmax = 850 mg / L (TDS), exceeding the Class III water standard (1000 mg / L) in an area of ​​0.45 km². B env = 1).

[0111] S3c Ecological Carrier Degradation Sub-model: Groundwater level drawdown: Theis formula was used. Q = 5000 m³ / d, T = 200 m² / d, S = 0.001. Wherein, Q This refers to the pumping flow rate. T The coefficient of conductivity is 1. S The storage coefficient is the maximum drawdown Δ over 5 years. h max = 45m, radius of influence R = 3500m.

[0112] Vegetation type: grassland (sheepgrass + needlegrass). z root = 0.8m, i fc = 0.28, E threshold = 0.55. The actual calculated maximum stress value. E stressmax = 0.72>0.55, the affected area is 12.5 km² ( B eco3D = 1).

[0113] S4. Boundary surface generation.

[0114] Comprehensive Boolean field: B forbid ( x,y,z,t ) = B geo ( x,y,z,t ) ∨B env ( x,y,z,t ) ∨Beco3D ( x,y,z,t The following area is designated as a no-mining zone: the area directly above and 200m in front of the working face. B geo = 1) Pollution protection zone: 0.45 km² downstream ( B env = 1) Ecologically confined areas: 12.5 km² of the Earth's surface ( B eco = 1) Safety boundary Ω safe ( t )for B forbid The complement of the set is used to extract isosurfaces via Marching Cubes. To address parameter uncertainties in multi-field coupling, Latin hypercube sampling is employed. V physchem and T process Generate a sample set of key parameters, run S3-S4 multiple times, count the frequency of each grid cell being marked as a no-entry zone, and extract the isosurface with 95% frequency as Ω. safe_95% ( t )boundary.

[0115] S5: Boundary output.

[0116] Generate a GeoTIFF layer and overlay it onto the mining area GIS platform.

[0117] The compliance report includes: Prohibited mining area: XYZ coordinate table (CGCS2000), area 0.32 km², validity period: permanent during mining; Pollution protection area: area 0.45 km², validity period: 10 years after mining; Ecological restriction area: area 12.5 km², validity period: from mining until water level recovery; Sensitivity: mining height ±0.5m results in H_f change ±12m, boundary displacement ±180m.

[0118] Example 2: Prediction of safety and environmental protection boundaries for shale gas development.

[0119] Project Background: A shale gas block is being developed using horizontal wells and staged fracturing. The horizontal stage is 1500m long and consists of 20 fracturing stages, with a single stage fluid volume of 1500m³ and a discharge rate of 12m³ / min. The target formation is 2500m deep and is overlain by a confined aquifer.

[0120] S1, Resource Form - Material Characteristics Encoding.

[0121] Resource type identifier: Shale gas development (primarily gas phase, with liquid phase fracturing). Retrieve physicochemical properties: Shale reservoir: r = 2450 kg / m³, s t= 4.5 MPa, s c = 65 MPa, gas compressibility factor Z = 1.12, gas diffusion coefficient D = 2.5 × 10⁻ 7 m² / s, Langmuir adsorption parameters L a = [ V L =3.2 m³ / t, P L =4.5 MPa], ion output coefficient I out = 0.35; fracturing fluid dynamic viscosity m = 5.2 mPa·s, saturation pressure P sat = 25 MPa, methane solubility S CH4 = 1200 mg / L, chloride ion concentration C Cl = 25000 mg / L (KCl solution), partial pressure of carbon dioxide P CO2 = 0.8 kPa.

[0122] Build V physchem = [2450, 4.5, 65, 1.12, 2.5e-7, 3.2, 4.5, 0.35, 5.2e-3, 25,1200, 25000, 0.8].

[0123] S2, Process parameter extraction.

[0124] T process ( t ) = [ P inj ( t ) =85MPa , Q inj ( t ) = 12 m³ / min, C frac ( t = [KCl, 2% gelling agent, 0.3% breaker], T reinj ( t =25℃, v des ( t ) = 2.5 × 10 4m³ / d, P vac ( t =50kPa, P back ( t =3.5MPa).

[0125] in, P inj ( t Inject pressure into the wellhead. Q inj ( t (This refers to the injection displacement.) C frac ( t This is a fracturing fluid formulation. T reinj ( t ( ) represents the temperature of the return fluid. v des ( t (This refers to the daily gas production.) P vac ( t () is the negative pressure of the vacuum pump. P back ( t () represents the back pressure at the wellhead.

[0126] S3, Multi-field Coupled Simulation.

[0127] S3a rock mass deformation-fracture sub-model: ABAQUS+Cohesive elements were used to simulate hydraulic fracture propagation. The fracture height was controlled by the stress difference between the upper and lower layers, with a simulated maximum fracture height of 85m, which did not penetrate the aquifer. B geo = 0).

[0128] S3b Groundwater Flow-Contaminant Transport Sub-Model: The fracturing fluid flowback rate is 30%, with the remaining 70% (approximately 21,000 m³) remaining in the formation. Migration of fracturing fluid additives (boron crosslinking agent, breaker residue) should be considered.

[0129] Simulation parameters: matrix permeability k m = 1×10⁻¹ 9 m², fracture permeability k f = 1×10⁻¹² m², crack porosity i f = 0.8.

[0130] A two-porosity model was adopted: matrix-fracture mass transfer coefficient α = 1×10⁻ 6 s⁻¹.

[0131] Results: At 10 years, the plume leading edge was 450m from the wellbore, and the Cl⁻ concentration was... C max = 1800 mg / L, exceeding the Class III water standard (250 mg / L) in an area of ​​0.18 km². B env = 1).

[0132] S3c Ecological Carrier Degradation Sub-model: Shale gas development has a relatively small impact on groundwater levels (Δh<2m), but surface fracturing operations occupy land and generate noise pollution.

[0133] Using a simplified model: E stress = f (Land occupation, noise, traffic) E stressmax = 0.25 < 0.55 ( B eco3D = 0).

[0134] S4. To address parameter uncertainties in multi-field coupling, Latin hypercube sampling is employed. V physchem and T process Generate a sample set of key parameters, run S3-S4 multiple times, count the frequency of each grid cell being marked as a no-entry zone, and extract the isosurface with 95% frequency as Ω. safe_95% ( t )boundary.

[0135] Boundary surface generation B forbid ( x,y,z,t )= B env ( B geo = B eco3D = 0) Safety boundary Ω safe ( t )for B env The supplementary set mainly constrains groundwater pollution protection.

[0136] S5. Boundary pollution protection range: area 0.18 km², validity period: 20 years after mining; key assumptions: fracture height is not out of control, interlayer integrity is maintained; sensitivity: ±10% of fracturing fluid volume results in ±15% pollution range.

[0137] Example 3: Prediction of Safety and Environmental Protection Boundaries for Hot Dry Rock Geothermal Development Project Background: A certain hot dry rock project adopts a dual-well system (one injection and one production), with a well depth of 4500m, a reservoir temperature of 180℃, an injection well flow rate of 80m³ / h, and a production well temperature drop limit of 20℃. The project has an operating life of 30 years.

[0138] S1, Resource Form - Material Characteristics Encoding.

[0139] Resource type identifier: Dry hot rock geothermal development (solid phase rock mass + liquid phase circulation).

[0140] Retrieve physicochemical properties: Granite reservoir: r = 2650 kg / m³, H = 6.5, s t = 8.5 MPa, σ c = 120MPa, thermal conductivity l = 2.8 W / (m·K), thermal diffusivity a = 1.2 × 10⁻ 6 m² / s, circulating working fluid (water): m = 0.18 mPa·s, r = 850 kg / m³, specific heat capacity c p = 4200 J / (kg·K).

[0141] V physchem = [2650, 6.5, 8.5, 120, 2.8, 1.2e-6, 0.18e-3, 850, 4200].

[0142] S2, Process parameter extraction.

[0143] T process ( t )= [ v circ ( t ) = 80 m³ / h, Δ T wp ( t =20℃, Y op ( t =30a).

[0144] in, v circ ( t ) represents the circulating flow rate, Δ T wp ( t ( ) represents the heat exchange temperature difference. Y op ( t() represents the system design life.

[0145] S3, Multi-field Coupled Simulation.

[0146] S3a rock mass deformation-fracture sub-model: Thermal extraction causes reservoir cooling and contraction, generating thermal stress.

[0147] A thermo-mechanical coupling model is used: s thermal = E·α T · Δ T / (1- n ); in, s thermal For thermal stress, E For elastic modulus, α T Δ is the coefficient of thermal expansion. T This refers to temperature changes.

[0148] Simulation results: After 30 years, the reservoir cooling radius is approximately 800m. Thermal stress leads to the propagation of microfractures and an increase in permeability of 2.3 times, but no large, interconnected fractures are formed. B geo = 0).

[0149] S3b Groundwater Flow-Contaminant Transport Sub-Model: The hot dry rock system is a closed loop with no external pollution sources. However, the following factors need to be considered: the water-rock reaction between the injected water and the surrounding rock, resulting in the dissolution of SiO2 and the release of trace heavy metals (As, F); and the impact of reinjection temperature changes on shallow groundwater.

[0150] Simulation results: The concentration of trace As increased from an initial 5 μg / L to 12 μg / L, which is lower than the Class III water standard (50 μg / L). B env = 0).

[0151] S3c Ecological Carrier Degradation Sub-model: The main impact of geothermal development on the surface ecosystem is land use (well sites, pipelines). A land use + noise model is used. E stressmax = 0.15 < 0.55 ( B eco3D = 0).

[0152] S4. To address parameter uncertainties in multi-field coupling, Latin hypercube sampling is employed. V physchem and T processA sample set of key parameters is generated, and S3-S4 are run multiple times to count the frequency at which each grid cell is marked as a forbidden zone. The isosurface with 95% frequency is extracted as the Ωsafe_95%(t) boundary. Boundary surface generation. B forbid = 0 (All areas are safe) B geo = B env = B eco3D = 0). Ω safe ( t The range is the entire reservoir, but the thermal decay zone (temperature drop > 10℃, radius 800m) needs to be marked.

[0153] S5. Boundary output thermal decay zone: radius 800m, validity period 30 years; Recommendation: after 30 years, the reservoir thermal recovery should be evaluated to determine whether to continue production; Sensitivity: injection flow rate ±10% results in thermal decay range ±8%.

[0154] It should be understood that the sequence number of each step in the above embodiments does not imply the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of the present invention.

[0155] The following are device embodiments of the present invention. For details not described in detail, please refer to the corresponding method embodiments described above.

[0156] Figure 2 A schematic diagram of the four-dimensional safety and environmental protection boundary prediction device for resource development provided in an embodiment of the present invention is shown. For ease of explanation, only the parts related to the embodiment of the present invention are shown, and are described in detail below: like Figure 2 As shown, the four-dimensional safety and environmental protection boundary prediction device 2 for resource development includes: The unified dimension conversion module 201 is used to obtain the materialized characteristic parameters corresponding to each resource type in the target development scheme and construct them into a unified materialized vector.

[0157] The process parameter acquisition module 202 is used to determine the process parameter tensor based on the sequence of mining process operation parameters of the target development plan.

[0158] The predicted spatiotemporal evolution acquisition module 203 is used to use the unified materialization vector and process parameter tensor as driving forces, and coupled the calculation of the spatiotemporal evolution results of equivalent plastic strain, pollutant concentration and vegetation stress index through the pre-established rock mass deformation-fracture sub-model, groundwater flow-pollutant transport sub-model and ecological carrier degradation sub-model.

[0159] The safe and environmentally friendly spatiotemporal boundary generation module 204 is used to generate a four-dimensional safe and environmentally friendly spatiotemporal boundary based on the spatiotemporal evolution results.

[0160] Figure 3 This is a schematic diagram of the four-dimensional safety and environmental protection boundary prediction device for resource development provided in an embodiment of the present invention. Figure 3 As shown, the four-dimensional safety and environmental protection boundary prediction device 3 for resource development in this embodiment includes a processor 300 and a memory 301. The memory 301 stores a computer program 302. When the processor 300 executes the computer program 302, it implements the steps in the above-described method embodiments. Alternatively, when the processor 300 executes the computer program 302, it implements the functions of each module / unit in the above-described device embodiments.

[0161] For the sake of simplicity and clarity, only the above-described functional modules / units are used as examples. In practical applications, the functions described above can be assigned to different functional modules / units as needed. These modules / units can be implemented in hardware, software, or a combination of both.

[0162] In the above embodiments, the descriptions of each embodiment have their own emphasis. Parts not detailed or described in a particular embodiment can be referred to in the relevant descriptions of other embodiments. Unless otherwise specified or in conflict with logic, the terminology and / or descriptions between different embodiments are consistent and can be referenced interchangeably. Technical features in different embodiments can be combined to form new embodiments based on their inherent logical relationships.

[0163] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention, and should all be included within the protection scope of the present invention.

Claims

1. A four-dimensional safety and environmental protection boundary prediction method for resource development, characterized in that, include: Obtain the materialized characteristic parameters corresponding to each resource type in the target development plan and construct them into a unified materialized vector; Based on the sequence of operational parameters for the mining process of the target development scheme, determine the process parameter tensor; Using the unified materialization vector and the process parameter tensor as driving forces, the spatiotemporal evolution results of equivalent plastic strain, pollutant concentration and vegetation stress index are coupled and calculated through the pre-established rock mass deformation-fracture sub-model, groundwater flow-pollutant transport sub-model and ecological carrier degradation sub-model. Based on the spatiotemporal evolution results, a four-dimensional safe and environmentally friendly spatiotemporal boundary is generated.

2. The four-dimensional safety and environmental protection boundary prediction method for resource development according to claim 1, characterized in that, The process uses the unified materialization vector and the process parameter tensor as driving forces, and couples the calculation of the spatiotemporal evolution results of equivalent plastic strain, pollutant concentration, and vegetation stress index through a pre-established rock mass deformation-fracture sub-model, groundwater flow-pollutant transport sub-model, and ecological carrier degradation sub-model, including: For each spatial discrete grid, the unified materialization vector and the process parameter tensor are input into the rock mass deformation-fracture sub-model, and the water-conducting fracture zone height, equivalent plastic strain and fracture permeability of the spatial discrete grid are output. The fracture permeability is updated to the groundwater flow-contaminant transport sub-model, and the unified materialization vector and the process parameter tensor are input into the updated groundwater flow-contaminant transport sub-model. The groundwater drawdown and rock pore pressure of the spatial discrete grid are output. When the height of the water-conducting fracture zone is greater than the fracture threshold, the contaminant concentration of the spatial discrete grid is output. The groundwater drawdown is updated to the ecological carrier degradation sub-model. The unified materialization vector and the process parameter tensor are input into the updated ecological carrier degradation sub-model, and the vegetation stress index of the spatial discrete grid is output. The rock mass pore pressure is updated to the rock mass deformation-fracture sub-model, and the iterative calculation for the next time point is initiated. Based on the time series values ​​of the equivalent plastic strain, pollutant concentration, and vegetation stress index corresponding to each spatial discrete grid, the spatiotemporal evolution results are determined.

3. The four-dimensional safety and environmental protection boundary prediction method for resource development according to claim 1, characterized in that, The generation of a four-dimensional safe and environmentally friendly spatiotemporal boundary based on the spatiotemporal evolution results includes: Based on the relationship between the equivalent plastic strain and the critical subsidence value within each spatial discrete grid, a time-varying geological hazard mining limit Boolean field is generated. Based on the relationship between the pollutant concentration in each spatial discrete grid and the preset pollution threshold, a groundwater environmental Boolean field that changes over time is generated. Based on the relationship between the vegetation stress index and the preset vegetation stress threshold in each spatial discrete grid, a three-dimensional ecological degradation dynamic Boolean field that changes over time is generated. Based on the geological hazard restricted mining Boolean field, the groundwater environmental protection Boolean field, and the three-dimensional ecological degradation restricted dynamic Boolean field, the four-dimensional safe and environmentally friendly spatiotemporal boundary is generated.

4. The four-dimensional safety and environmental protection boundary prediction method for resource development according to claim 3, characterized in that, The generation of the four-dimensional safe and environmentally friendly spatiotemporal boundary based on the geological hazard-limited mining Boolean field, the groundwater environmental protection Boolean field, and the three-dimensional ecological degradation-limited dynamic Boolean field includes: By performing a logical OR operation on the geological hazard restricted mining Boolean field, the groundwater environmental protection Boolean field, and the three-dimensional ecological degradation restricted movement Boolean field, a comprehensive prohibited entry Boolean field is obtained; The four-dimensional safe and environmentally friendly spacetime boundary is determined based on the complement of the comprehensive forbidden Boolean field.

5. The four-dimensional safety and environmental protection boundary prediction method for resource development according to any one of claims 1-4, characterized in that, The formulas for the rock mass deformation-fracture sub-model include: ∇ ·σ + ρ g b = 0 ； σ = ( 1 - D )· C : ε e ； Ḋ = f ( ε eq , σ 1 , σ 3 , M, H f )； Hf =∑ n i=1 hi ⋅ I ( Si < Scr )； k f = k 0 · exp ( α · ε v ) ( ε v > ε cr )； in, σ For rock mass stress tensor, ρ g The rock mass density is driven by the unified materialization vector. b For the volume forces of the rock mass, ∇ · For divergence operators, D For rock mass damage variables, C Let be the rock mass elastic stiffness tensor. ε e Let be the elastic strain tensor of the rock mass, and be the double dot product. ε eq For the equivalent plastic strain of the rock mass, σ 1 The maximum principal stress of the rock mass, σ 3 The minimum principal stress of the rock mass, M The sampling height is driven by the process parameter tensor. H f The height of the water-conducting fracture zone, hi For the first i Thickness of rock layers, Si For the first i The drop in rock strata value, Scr This represents the critical subsidence value of the rock mass. I For indicator functions, n This represents the total number of overburden layers. k f The fracture permeability, k 0 The initial permeability of the rock mass. α The damage-permeability coupling coefficient is... ε v For the volumetric strain of the rock mass, exp ( · ) is an exponential function, ε cr The critical volumetric strain of the rock mass is given.

6. The four-dimensional safety and environmental protection boundary prediction method for resource development according to any one of claims 1-4, characterized in that, The formulas for the groundwater flow-contaminant transport sub-model include: ∂ ( θp ) / ∂t = ∇⋅( k f / μ⋅ ∇ p ) +Q s ； v = -k f / μ· (∇ p + ρ y g ∇ z )； ∂ ( θC ) / ∂t = ∇·( θD· ∇ C ) - ∇( θvC ) + θR ( C ) + q s C s ； in, θ The rock mass porosity is driven by the unified materialization vector. p The pore pressure of the rock mass is... k f The fracture permeability, μ For the unified materialization vector-driven hydrodynamic viscosity, Q s Let ∇ be the source and sink terms of the seepage equation driven by the unified materialization vector, ∇ be the gradient operator, and ∇⋅ be the divergence operator. v The average flow velocity of the fluid in the void. ρ y For the fluid density driven by the unified materialization vector, g It is the acceleration due to gravity. z For elevation, C The concentration of the pollutant, D The hydrodynamic dispersion coefficient tensor driven by the unified materialization vector. R ( C () represents a chemical reaction term. C s The source and sink concentrations are driven by the process parameter tensor. q s These are the source and sink terms of the concentration equation driven by the unified materialization vector.

7. The four-dimensional safety and environmental protection boundary prediction method for resource development according to any one of claims 1-4, characterized in that, The formula for the ecological carrier degradation sub-model includes: E stress ( x,y,t ) =w 1 ·f 1 (D h / h max ) + w 2 ·f 2 (D θ r / θ fc ) +w 3 ·f 3 ( PET / P ); in, E stress ( x,y,t ) for time t Plane coordinates ( x , y The vegetation stress index described below, Δh For the aforementioned groundwater drawdown, h max The maximum root water absorption depth of vegetation driven by the unified materialization vector. θ r For the residual content, θ fc For the unified materialized vector-driven field water holding capacity, PET For potential evaporation, P For precipitation, f 1 The first vegetation type response coefficient is driven by the unified materialization vector. f 2 The second vegetation type response coefficient is driven by the unified materialization vector. f 3 The response coefficient of the third vegetation type driven by the unified materialization vector. w 1 The first weight coefficient driven by the unified materialization vector. w 2 The second weighting coefficient is driven by the unified materialization vector. w 3 The third weight coefficient driven by the unified materialization vector satisfies w 1 + w 2 + w 3 =1.

8. The four-dimensional safety and environmental protection boundary prediction method for resource development according to any one of claims 1-4, characterized in that, The process of obtaining the materialized characteristic parameters corresponding to each resource type in the target development plan and constructing them into a unified materialized vector includes: Obtain the physicochemical property parameters corresponding to the solid phase resources, liquid phase resources and gas phase resources in the target development scheme; The physicochemical property parameters corresponding to the solid phase resources, liquid phase resources and gas phase resources are mapped to a unified physical property coordinate system for standardization to obtain standardized solid phase parameters, standardized liquid phase parameters and standardized gas phase parameters. The standardized solid-phase parameters, the standardized liquid-phase parameters, and the standardized gas-phase parameters are concatenated and multiplied by the transformation matrix to obtain the unified physicalization vector.

9. A four-dimensional safety and environmental protection boundary prediction device for resource development, characterized in that, include: The unified dimension conversion module is used to obtain the materialized characteristic parameters corresponding to each resource type in the target development plan and construct them into a unified materialized vector. The process parameter acquisition module is used to determine the process parameter tensor based on the sequence of mining process operation parameters of the target development scheme; The predicted spatiotemporal evolution acquisition module is used to use the unified materialization vector and the process parameter tensor as driving forces, and through the pre-established rock mass deformation-fracture sub-model, groundwater flow-pollutant transport sub-model and ecological carrier degradation sub-model, to couple and calculate the spatiotemporal evolution results of equivalent plastic strain, pollutant concentration and vegetation stress index. The safe and environmentally friendly spatiotemporal boundary generation module is used to generate a four-dimensional safe and environmentally friendly spatiotemporal boundary based on the spatiotemporal evolution results.

10. A four-dimensional safety and environmental protection boundary prediction device for resource development, characterized in that, It includes a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the four-dimensional safety and environmental protection boundary prediction method for resource development as described in any one of claims 1 to 8.

Images