Coalfield aquifer hydrological parameter optimization and numerical simulation method
By combining sequential instruction simulation and physical information neural network optimization methods with a distributed optical fiber sensing system, the problems of accuracy and dynamic change in obtaining hydrological parameters of coalfield aquifers were solved, achieving high-precision simulation of groundwater flow and identification of water control target areas.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- GEOPHYSICAL SURVEY TEAM OF SHANDONG COALFIELD GEOLOGY BUREAU
- Filing Date
- 2026-05-13
- Publication Date
- 2026-07-24
AI Technical Summary
Existing technologies struggle to accurately obtain hydrological parameters of coalfield aquifers, especially the spatial distribution and dynamic changes of permeability coefficient, water storage coefficient, and porosity. This leads to a high risk of mine water inrush during coal mining. Traditional numerical simulation methods are costly and have low accuracy, failing to reflect the spatial heterogeneity and dynamic changes of these parameters.
By employing sequential indicator simulation and physical information neural network global parameter optimization methods, combined with real-time monitoring of temperature and strain data by a distributed fiber optic sensing system, and through Bayesian inversion and adaptive weight adjustment, hydrological parameters are optimized to construct a high-precision groundwater flow model, thereby realizing the inversion and simulation of dynamic parameters.
It achieves spatial continuity, temporal dynamics, and high-resolution inversion of hydrological parameters of coalfield aquifers, solves the problem of inaccurate parameters in traditional models, improves simulation accuracy and parameter optimization efficiency, and meets the real-time flow field prediction needs of coal mine water control.
Smart Images

Figure CN122287383B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of coalfield hydrogeological exploration and groundwater numerical simulation technology, and more specifically, to a method for optimizing and numerically simulating hydrological parameters of coalfield aquifers. Background Technology
[0002] Coal is my country's primary energy source, and safe coal mine production is a crucial guarantee for national energy security. During coal mining, mine water inrush is one of the major hazards threatening safe production. The root cause lies in a lack of understanding of the hydrogeological conditions of coalfield aquifers, particularly insufficient knowledge of the spatial distribution and dynamic changes of key hydrological parameters (such as permeability coefficient, water storage coefficient, and porosity).
[0003] Currently, the main techniques for obtaining and numerically simulating hydrological parameters of coalfield aquifers are as follows: hydrological parameters are derived by combining single-hole or multi-hole pumping tests with analytical methods. This method is costly and time-consuming, and the obtained parameters only reflect the average characteristics within a small area near the borehole, failing to characterize the spatial heterogeneity of the parameters. Laboratory experiments, due to sample disturbance and scale effects, cannot accurately reflect the seepage characteristics of aquifers under in-situ conditions. Summary of the Invention
[0004] To overcome the aforementioned deficiencies of existing technologies, this invention provides a method for optimizing and numerically simulating hydrological parameters of coalfield aquifers. Through sequential indicator simulation and global parameter optimization using physical information neural networks, it fully characterizes the strong heterogeneity and anisotropy of coalfield aquifers, fundamentally solving the fatal flaws of inaccurate parameters and disconnection from geological reality in traditional numerical models.
[0005] To achieve the above objectives, the present invention provides the following technical solution: A method for optimizing and numerically simulating hydrological parameters of a coalfield aquifer includes the following steps: Step S1: Acquire drilling monitoring data and historical hydrogeological data from multiple boreholes within the target coalfield area, construct an initial aquifer structure model, and generate an initial hydrological parameter field based on the initial aquifer structure model; Step S2: Based on the initial hydrological parameter field, construct an initial groundwater numerical simulation model, and calibrate and verify the initial groundwater numerical simulation model using preset observation well water level data to obtain a baseline numerical simulation model; Step S3: Deploy a distributed fiber optic sensing system in the target coalfield area to acquire temperature data at multiple spatial locations in real time. Step S4: Construct a parameter optimizer based on a physical information neural network. The optimized parameter set is generated by taking the simulated water level of the benchmark numerical simulation model, the dynamic hydrological parameter field, and the water level data from the observation wells as inputs. The optimized parameter set is generated by iteratively optimizing the hydrological parameters by solving the loss function embedded in the partial differential equation of groundwater flow. Step S5: Input the optimized hydrological parameter set into the benchmark numerical simulation model, update the model parameters, and perform multi-scale simulation calculations to output high-precision prediction results of aquifer water level and flow field distribution.
[0006] In a preferred embodiment, the construction of the initial aquifer structure model in step S1 specifically involves: cleaning and standardizing the drilling monitoring data and historical hydrogeological data to extract lithological stratification information and marker layer depths for each borehole; using a sequential indicator simulation algorithm, with the lithological stratification information as hard data and the stratigraphic data interpreted from seismic exploration as soft data, performing three-dimensional spatial interpolation and random simulation to generate multiple equally probable lithological distribution realizations; calculating the average probability distribution of the multiple equally probable lithological distribution realizations and using it as the initial aquifer structure model.
[0007] In a preferred embodiment, the calibration and verification step S2 specifically involves: extracting the first time period data from the observation well water level data as calibration period data and the second time period data as verification period data; using a parallel particle swarm optimization algorithm to automatically calibrate the hydrological parameters of the initial groundwater numerical simulation model with the goal of minimizing the root mean square error between the simulated water level and the measured water level in the calibration period data; using the verification period data to verify the simulation accuracy of the calibrated model; and determining the current model as the benchmark numerical simulation model when the Nash efficiency coefficient in the verification period is greater than a preset threshold.
[0008] In a preferred embodiment, the inversion of the dynamic hydrological parameter field in step S3 specifically involves: preprocessing the temperature and strain time series data by denoising and normalization to extract its spatiotemporal variation characteristics; constructing an inversion model based on the heat conduction equation and the thermo-solid coupling constitutive relation, using the preprocessed temperature and strain time series data as input to the inversion model, and using the aquifer permeability coefficient and porosity as the parameters to be inverted; using a Bayesian inversion framework and combining it with the Markov chain Monte Carlo method to solve the inversion model and obtain the posterior probability distribution of the parameters to be inverted; and extracting the mean of the posterior probability distribution as the dynamic hydrological parameter field.
[0009] In a preferred embodiment, step S4, which involves solving the loss function embedded in the partial differential equation of groundwater flow, specifically involves: constructing a deep neural network consisting of an input layer, multiple fully connected hidden layers, and an output layer. The input layer receives spatial and temporal coordinates, and the output layer outputs the predicted hydraulic head value. A total loss function is constructed, including data fitting terms, initial condition constraints, boundary condition constraints, and physical equation constraints. The physical equation constraints are based on the groundwater flow continuity equation and Darcy's law, and are used to constrain the predicted hydraulic head value output by the deep neural network to satisfy the physical laws of groundwater flow. An adaptive moment estimation optimizer is used to iteratively update the network weights of the deep neural network with the goal of minimizing the total loss function. After training, optimized hydrological parameters are extracted from the deep neural network.
[0010] In a preferred embodiment, the multi-scale simulation calculation in step S5 specifically involves: constructing a nested grid system, using locally refined grids in the well periphery and geologically complex areas, and coarse grids in other areas; mapping the optimized hydrological parameter set onto each grid cell of the nested grid system; discretizing the groundwater flow control equations using the finite volume method based on the nested grid system to construct a large-scale sparse linear equation set; and solving the large-scale sparse linear equation set in parallel using the algebraic multigrid preprocessing conjugate gradient method to obtain high-precision prediction results of aquifer water level and flow field distribution.
[0011] In a preferred embodiment, the method further includes step S6: performing three-dimensional visualization rendering of the high-precision aquifer water level and flow field distribution prediction results to generate a dynamic interactive display interface; identifying and delineating water-rich anomaly areas and potential water-conducting channels based on the dynamic interactive display interface; and outputting the boundary information of the water-rich anomaly areas and potential water-conducting channels as the target area coordinate set for coal mine water control work.
[0012] In a preferred embodiment, the identification and delineation of water-rich anomaly zones and potential water-conducting channels specifically involves: calculating the spatial variability of the hydraulic gradient and permeability coefficient tensors in the high-precision aquifer level and flow field distribution prediction results; calculating the flow dominance factor for each grid cell based on the hydraulic gradient and permeability coefficient tensors, whereby the flow dominance factor characterizes the likelihood of preferential water flow; performing three-dimensional spatial clustering analysis on the flow dominance factor to identify continuous regions with high flow dominance factors as potential water-conducting channels; and identifying regions with hydraulic gradients below a preset threshold as water-rich anomaly zones.
[0013] In a preferred embodiment, the parameter optimizer based on the physical information neural network further includes a step of dynamically adjusting the weights of the physical equation constraint terms in the loss function during the iterative optimization process. Specifically, this involves: calculating the average absolute error between the simulated water level and the measured water level in the current iteration step; adaptively adjusting the weight coefficients of the physical equation constraint terms according to the changing trend of the average absolute error; and increasing the weight coefficients of the physical equation constraint terms when the average absolute error tends to stabilize, so as to strengthen the constraints of physical laws.
[0014] The technical effects and advantages of the method for optimizing and numerically simulating hydrological parameters of coalfield aquifers in this invention are as follows: This invention overcomes the limitations of traditional pumping / pressure water tests, which suffer from single-point discreteness, time lag, and the inability to obtain only static averages, by using distributed fiber optic temperature and strain dual-field monitoring and coupled Bayesian inversion of heat, solid, and fluid. It achieves spatial continuity, temporal dynamics, and high-resolution inversion of core hydrological parameters such as aquifer permeability and porosity, accurately capturing the dynamic process of rock mass fracture development and parameter mutations under coal mining disturbances. At the same time, through sequential indicator simulation and physical information neural network (PINN) global parameter optimization, it fully characterizes the strong heterogeneity and anisotropy of coalfield aquifers, fundamentally solving the fatal defects of inaccurate parameters and disconnection from geological reality in traditional numerical models.
[0015] This invention embeds the groundwater flow control equations into PINN and employs an adaptive weight adjustment strategy, completely resolving the dual industry dilemmas of traditional pure data-driven models (black box nature, lack of physical meaning, susceptibility to overfitting, poor generalization) and pure physics-driven models (inability to adapt to dynamic changes and susceptibility to local optima). It ensures the dynamic accuracy of parameters through real-time monitoring data and implements global constraints through hydrogeological and physical laws, upgrading parameter optimization from manual trial-and-error fitting to global intelligent optimization, improving parameter optimization efficiency by more than an order of magnitude, while simultaneously guaranteeing the physical rationality and engineering credibility of the optimization results.
[0016] By employing a proprietary scheme involving nested local mesh refinement, finite volume method conservation discretization, and parallel solution of conjugate gradients through algebraic multigrid preprocessing, ultra-high precision characterization of core water control areas such as well perimeter and geologically complex zones is achieved while controlling the overall computational load. Simultaneously, the simulation time for large-scale coalfield models is reduced from several days to minutes, completely resolving the deadlock of insufficient precision in traditional uniform mesh with coarse mesh and explosive computational load in fine mesh, thus meeting the real-time flow field prediction requirements under dynamic changes in coal mining.
[0017] Through a complete process including full-element 3D dynamic interactive visualization, original quantitative indicators of water flow advantage factors, 3D spatial clustering risk identification, and standardized output of coordinate sets for water control target areas, discrete numerical results that were originally only interpretable by hydrogeological professionals are transformed into engineering results that can be understood by all positions in coal mines and directly applied to underground construction. This breaks down professional barriers and achieves a seamless connection between high-precision simulation, accurate risk identification, and targeted engineering treatment, allowing numerical simulation results to truly serve as a basis for on-site decision-making in coal mine water control.
[0018] From the realization of multiple equal probability lithological distributions in the initial modeling stage, to the solution of Bayesian posterior probability distributions in the parameter inversion stage, to the global constraint of physical laws in the parameter optimization stage, the inherent uncertainty of coalfield hydrogeological modeling and simulation is quantified and controlled throughout the entire process. This avoids the one-sidedness and subjective bias of traditional single deterministic models, and fundamentally improves the reliability and reproducibility of simulation results, fully meeting the high reliability requirements of coal mine safety production for decision-making basis. Attached Figure Description
[0019] Figure 1 The overall flowchart of the method for optimizing hydrological parameters of coalfield aquifers and numerical simulation provided in the embodiments of the present invention is shown.
[0020] Figure 2 This is a schematic diagram of a nested grid system in an embodiment of the present invention. Detailed Implementation
[0021] 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 of ordinary skill in the art without creative effort are within the scope of protection of the present invention.
[0022] Example 1, Figure 1 and Figure 2This invention presents a method for optimizing and numerically simulating hydrological parameters of a coalfield aquifers, comprising the following steps: Step S1: Acquire drilling monitoring data and historical hydrogeological data from multiple boreholes within the target coalfield area, construct an initial aquifer structure model, and generate an initial hydrological parameter field based on the initial aquifer structure model; Step S2: Based on the initial hydrological parameter field, construct an initial groundwater numerical simulation model, and calibrate and verify the initial groundwater numerical simulation model using preset observation well water level data to obtain a baseline numerical simulation model; Step S3: Deploy a distributed fiber optic sensing system in the target coalfield area to acquire multiple spatial parameters in real time. The temperature and strain time series data of the location are used to invert the dynamic hydrological parameter field based on the temperature and strain time series data; Step S4: Construct a parameter optimizer based on physical information neural network, take the simulated water level, dynamic hydrological parameter field and observation well water level data of the benchmark numerical simulation model as input, solve the loss function embedded in the partial differential equation of groundwater flow, iteratively optimize the hydrological parameters, and generate the optimized hydrological parameter set; Step S5: Input the optimized hydrological parameter set into the benchmark numerical simulation model, update the model parameters, and perform multi-scale simulation calculations to output high-precision prediction results of aquifer water level and flow field distribution.
[0023] In this embodiment, the initial aquifer structure model is constructed in step S1, specifically by: cleaning and standardizing the drilling monitoring data and historical hydrogeological data, extracting the lithological stratification information and marker layer depth of each borehole; using a sequential indicator simulation algorithm, with lithological stratification information as hard data and seismic exploration interpretation stratigraphic data as soft data, performing three-dimensional spatial interpolation and random simulation to generate multiple equally probable lithological distribution realizations; calculating the average probability distribution of the multiple equally probable lithological distribution realizations, and using it as the initial aquifer structure model.
[0024] It should be noted that the two types of data sources processed in this step are the two most complementary and readily available core geological data in a coalfield scenario, and also the foundational anchor points for modeling: Drilling monitoring data: This refers to real-time logging data (gamma, resistivity, density, sonic transit time, etc.), well logging data (cuttings logging, drilling time logging, gas logging), and core logging data acquired during drilling operations within the target coalfield area. The core characteristics of this type of data are extremely high vertical resolution, strong real-time performance, and high lithological identification accuracy, serving as the core basis for vertical stratigraphic division in single boreholes. Historical hydrogeological data: This refers to existing regional geological survey reports, complete data from previous exploration boreholes, hydrogeological test results, stratigraphic correlation reports, regional structural maps, etc., for the target area. The core characteristics of this type of data are wide coverage and mature understanding of regional geological patterns, which can compensate for the insufficient number of newly drilled boreholes.
[0025] The core steps in resolving the issues of chaotic, fragmented, erroneous, and incompatible multi-source geological data in coalfields, and also the most easily overlooked crucial steps in traditional modeling, are: Data Cleaning: The core is to remove invalid, erroneous, and abnormal data, and to supplement key missing information. Specifically, this includes: removing logging data jumps caused by stuck drill bits, drill string drops, or instrument malfunctions during drilling; correcting abnormal data in historical data that are manually recorded errors or have mismatched coordinates / depths; supplementing missing marker layer depth data for some boreholes using regional stratigraphic correlation; and standardizing the borehole elevation and vertical depth benchmarks for all boreholes to eliminate depth conversion errors for inclined and directional boreholes. Standardization Processing: The core is to resolve the issue of inconsistent data calibers across multiple sources, which is a prerequisite for subsequent 3D spatial modeling. The most crucial aspect is the standardization of the lithological classification system: the lithological cataloging standards vary greatly across different periods and units in coalfields (e.g., some data merge fine sandstone and siltstone into a single name, while others are finely separated). This step requires establishing a unified lithological coding system adapted to the hydrogeological scenarios of coalfields, with the core classification into two main categories: aquifer lithology (coarse sandstone, medium sandstone, conglomerate, and other highly permeable layers) and relatively impermeable lithology (mudstone, silty mudstone, coal seams, and other weakly permeable / impermeable layers). This involves unifying the lithological naming rules and grain size classification standards; simultaneously, unifying the spatial coordinate system and elevation datum of all data to completely eliminate spatial misalignment issues.
[0026] Lithological stratification information extraction: Accurate extraction of lithological interface depth, single-layer lithology type, single-layer thickness, and stratigraphic contact relationships along the vertical direction for each borehole. This information is the hard data (measured deterministic data) for subsequent modeling, serving as the measured anchor points for the 3D model and ensuring that the geological attributes of the model at the borehole location are 100% consistent with reality. Marker layer depth extraction: Marker layers are strata with stable distribution within the coalfield area, easily identifiable lithological characteristics, minimal thickness variation, and strong regional comparability (such as the roof and floor of the main coal seam, bauxite mudstone, oolitic limestone, etc.). They serve as a benchmark for regional stratigraphic correlation. Extracting marker layer depth can resolve stratigraphic misalignment issues caused by tectonic movements such as folding and faulting between different boreholes, ensuring that the stratigraphic framework of the 3D model conforms to regional tectonic patterns and fundamentally avoiding misalignment of aquifer roof and floor boundaries due to stratigraphic correlation errors.
[0027] Traditional coalfield modeling often uses raw borehole data directly, ignoring the inconsistencies in coordinate systems and lithological standards across multiple data sources. This leads to misaligned stratigraphic frameworks and chaotic lithological classifications, resulting in initial models that fundamentally do not reflect geological reality. No matter how hydrological parameters are subsequently adjusted, reliable simulation results cannot be obtained. This step transforms the chaotic, heterogeneous multi-source data into a unified, reliable, and standardized dataset suitable for 3D modeling, laying a solid data foundation for subsequent modeling.
[0028] Traditional coalfield modeling commonly uses Kriging interpolation, but it has an intractable fatal flaw: Kriging is an interpolation method for continuous variables, while lithology is a typical discrete categorical variable (sandstone and mudstone are categorical attributes, not continuous values). Using Kriging to interpolate lithology will result in serious distortion, and the interpolation results will be overly smoothed, making it completely unable to characterize the strong heterogeneity of coalfield strata (such as sandstone lenses and lithological facies transition zones, which are the core hidden dangers of coal mine water inrush).
[0029] Sequential Indicator Simulation (SIS) is a stochastic simulation method developed specifically for discrete categorical variables in the field of geostatistics. It is perfectly suited for coalfield lithology modeling scenarios, and its core advantages include: complete adaptation to the discrete categorical attributes of lithology, with interpolation results conforming to geological sedimentary laws; it not only strictly adheres to the hard data measured in boreholes, but also restores the discreteness and heterogeneity of lithology spatial distribution without over-smoothing, and can accurately characterize key geological bodies such as aquifer lenses and phase transition zones; it can generate multiple lithology distribution models with equal probability, rather than a single deterministic model, realizing the quantitative characterization of uncertainty in geological modeling and solving the industry problem that traditional single models cannot assess the reliability of results.
[0030] While borehole data offers high precision, it is sparse (often spaced hundreds or even thousands of meters apart), making it impossible to control lithological variations between boreholes using only borehole data. Seismic exploration data covers the entire area and is highly continuous, but it represents indirect interpretation and is subject to some ambiguity. This step perfectly resolves this contradiction by fusing hard and soft data: Hard data refers to the borehole lithological stratification information extracted in the first step, which is 100% measured and absolutely reliable point data. In sequential indicator simulations, the lithological values at the hard data locations remain constant, and the simulation results must strictly adhere to the hard data to ensure the absolute accuracy of the model at the borehole locations. Soft data refers to the stratigraphic data interpreted from seismic exploration, which is regionally continuous, constrains stratigraphic trends, and is not directly measured surface / volume data. Seismic exploration can achieve full coverage of the target area and accurately interpret the stratigraphic boundaries, undulations, and spatial trends of lithological facies changes, effectively filling information gaps between boreholes. However, because it is an indirect interpretation, it has a certain degree of ambiguity. Therefore, as soft data, it is used to constrain the lithological distribution trends between boreholes and avoid excessive randomness in simulation results.
[0031] Geostatistical Feature Analysis: A variogram analysis is performed on the standardized lithological data to quantify the spatial correlation characteristics of different lithologies. Key parameters of the variogram (range, sill value, nugget value) are calculated. In simpler terms, this clarifies the spatial range within which a certain type of lithology exhibits correlation (e.g., the correlation range for sandstone along the bedding plane is 500m, and vertically it is 50m, meaning that lithology within 500m of the bedding plane is likely to remain consistent; beyond this distance, the correlation decreases significantly). This forms the statistical basis for stochastic simulations. 3D Simulation Grid Construction: Based on the planar extent and vertical stratigraphic thickness of the target coalfield, a 3D grid adapted for hydrological simulation is constructed (e.g., a 50m×50m planar grid and a 2m vertical grid). This balances simulation accuracy and computational efficiency. The grid cell is the smallest unit for subsequent lithological assignment, hydrological parameter assignment, and numerical simulation.
[0032] Sequential random simulation traversal: Each node of the 3D mesh is traversed through a random path: If the node has borehole hard data, the corresponding lithology is directly assigned and kept fixed in subsequent simulations; if the node does not have hard data, the conditional probability of the node belonging to different lithology types is calculated based on the surrounding nodes that have already been assigned values (hard data + nodes that have been simulated), the pre-calculated variogram, and the trend constraints of the seismic horizon soft data; through random sampling, the lithology of the node is assigned according to the calculated conditional probability.
[0033] Multiple realizations are generated: Each time a full grid traversal simulation is completed, one lithological distribution realization with equal probability is generated (i.e., one three-dimensional lithological distribution model that conforms to the hard data of boreholes, the soft data of seismic earthquakes, and the laws of geological statistics); by repeating the simulation multiple times (usually 50-100 times), dozens of lithological distribution realizations with equal probability can be generated.
[0034] This method completely solves the problems of interpolation distortion and over-smoothing of discrete lithology variables in traditional interpolation methods, and accurately characterizes the strong heterogeneity of coalfield aquifers. By integrating hard and soft data, it takes into account both the absolute accuracy of borehole points and the overall trend of regional strata, and solves the problem of unconstrained lithology distribution among boreholes caused by the sparseness of coalfield boreholes. Through multiple equally probable generation, it quantifies the inherent geological uncertainty of aquifer structure modeling for the first time, and provides a foundation for the reliability assessment of subsequent simulation results.
[0035] For each grid cell in the 3D mesh, the frequency of occurrence of a certain lithology in all equally probable lithology distribution realizations is statistically analyzed, and this frequency is used as the probability value of the corresponding lithology for that cell. For example, if 100 equally probable realizations are generated, and a certain grid cell is assigned the value of aquifer sandstone in 85 realizations and mudstone in 15 realizations, then the probability of sandstone in that cell is 85%, and the probability of mudstone is 15%. This process is repeated for all cells in the entire 3D mesh, ultimately yielding the average probability distribution model of various lithologies in the 3D space of the entire study area.
[0036] While single-lithological distributions conform to existing data constraints, they exhibit strong randomness and cannot represent the lithological distributions that best reflect geological realities. In contrast, the average probability distribution integrates geological information from all randomly generated distributions, minimizing the random bias of single simulations. It is the most representative and geologically consistent result among all equally probable distributions. The average probability distribution is not a traditional, black-and-white deterministic lithological classification; rather, it preserves the lithological probability attribute of each grid cell, quantifying the reliability of whether a location is an aquifer / impermeable layer. This provides continuous, geologically consistent constraints for subsequent hydrological parameter assignments.
[0037] Based on the average probability distribution, the spatial division of the aquifer and impermeable layer is completed, and the final three-dimensional initial aquifer structure model is constructed: a lithological probability threshold is set, and combined with coalfield hydrogeological experience, aquifer and impermeable layer units are divided: for example, grid units with a sandstone probability ≥60% are divided into aquifer units, and grid units with a mudstone probability ≥80% are divided into impermeable layer units; combined with the marker layer depth constraint extracted in the first step, the top and bottom plate boundaries, spatial distribution range, thickness variation, and spatial connectivity of the aquifer are accurately determined; finally, a three-dimensional aquifer structure model with lithological probability attributes is output. This model not only clarifies the spatial morphology of the aquifer, but also quantifies the lithological reliability of each location, serving as the geological carrier for all subsequent hydrological simulation work.
[0038] This approach balances the characterization of geological heterogeneity in stochastic simulation with the reliability of the results, avoiding the limitations of a single model. The output structural model with probabilistic attributes provides direct geological constraints for the generation of the initial hydrological parameter field in step S1: units with higher sandstone probability are assigned larger values for hydrological parameters such as permeability and porosity, while units with higher mudstone probability are assigned smaller values. This ensures that the initial hydrological parameter field fully conforms to the spatial distribution of lithology from the outset, completely avoiding the irrationality of layered homogeneous assignments in traditional modeling. It fundamentally improves the rationality of the initial parameter field, reduces the number of iterations for subsequent model calibration and parameter optimization, and significantly improves the convergence speed and accuracy of the final simulation.
[0039] In this embodiment, calibration and verification are performed in step S2, specifically as follows: the first time period data from the observation well water level data is extracted as the calibration period data, and the second time period data is used as the verification period data; a parallel particle swarm optimization algorithm is used to automatically calibrate the hydrological parameters of the initial groundwater numerical simulation model with the goal of minimizing the root mean square error between the simulated water level and the measured water level in the calibration period data; the simulation accuracy of the calibrated model is verified using the verification period data, and when the Nash efficiency coefficient in the verification period is greater than a preset threshold, the current model is determined as the benchmark numerical simulation model.
[0040] It should be noted that the water level data used in the observation wells is the only continuous time-series measured true value that can directly reflect the dynamic changes of groundwater in the coalfield scenario. It is fundamentally different from the drilling and historical geological data mentioned above. It comes from long-term hydrological observation wells specially deployed in the coalfield. It is a continuous water level monitoring sequence on a minute / hour / day scale. It fully records the dynamic changes of groundwater under all working conditions, including coal mining, atmospheric precipitation infiltration, aquifer dewatering, and overflow recharge. It is the only and irreplaceable measured constraint basis for model calibration.
[0041] This step addresses the fundamental problem of overfitting, a long-standing and fatal flaw in coalfield hydrological models. Traditional modeling often uses all water level data for fitting, resulting in a model that can only reproduce known data and becomes completely distorted when the time period changes, lacking the predictive capabilities necessary for water control. This step avoids overfitting risks from the data source by strictly defining independent and non-overlapping time series. Its coalfield-specific segmentation rules are as follows: Calibration Period (First Time Period): Prioritizes long-term time series data with complete hydrological dynamic characteristics, good data continuity, and no extreme anomalies. It is mandatory to cover a complete hydrological year (including the wet season, normal season, and dry season), and must include typical working conditions such as normal coal mine drainage, working face advancement, and seasonal shutdown, with a duration of no less than 70% of the total observation sequence. Its core function is to provide sufficiently rich dynamic constraints for model parameter calibration, allowing the model to fully learn the core laws of groundwater flow in the study area, rather than just fitting water level values at a single point or under a single working condition. The validation period (second time period) must be a continuous period completely independent of the rate period, with no data overlap. Priority should be given to periods following the rate period, with a duration of no less than 30% of the total observation sequence. It is mandatory to include hydrological conditions that differ from the rate period (e.g., the rate period is normal mining, while the validation period includes extreme conditions such as water inrush at the working face, large-scale mining shutdown, and concentrated rainstorm infiltration). Its core function is to test the model's generalization and predictive ability, not to reproduce known data. It verifies whether the model can still accurately predict groundwater dynamics under conditions not included in calibration. Dedicated preprocessing is provided: For well water level data, abnormal jump values caused by temporary forced drainage, water inrush accidents, and human monitoring failures are specifically removed. Short-term data gaps of no more than 3 days are filled in using linear interpolation. The time step and elevation benchmark are unified between the rate period and the validation period to ensure the consistency and comparability of the time series data.
[0042] Traditional coalfield hydrological models generally lack rigorous independent validation processes, relying solely on visual fitting to judge model effectiveness. This results in many models performing well in laboratory settings but making completely inaccurate predictions in the field, failing to provide a reliable basis for coal mine water control. This step, through standardized independent time-series partitioning, establishes an isolation mechanism between model fitting and prediction at the data level, completely avoiding the risk of overfitting and ensuring the predictive reliability of the final benchmark model.
[0043] The parameters to be calibrated in this step are the core sensitive hydrological parameters that determine the accuracy of groundwater simulation. These include the horizontal / vertical permeability coefficient, porosity, specific yield / storage coefficient, and overflow coefficient of each aquifer. These parameters are characterized by strong spatial heterogeneity, high dimensionality, and nonlinear coupling.
[0044] Traditional coalfield model calibration generally adopts a manual trial-and-error method, which relies entirely on the experience of hydrogeologists. After manually adjusting parameters, the model is repeatedly run to check the fitting effect. This method has fatal flaws that cannot be solved: it can only adjust 3-5 core parameters and cannot achieve global optimization of multiple parameters; it is extremely inefficient, with the calibration cycle of a medium-sized coalfield model taking several months; it is very easy to get trapped in local optima; and the results of calibration by different personnel vary greatly, with no unified quantitative standard.
[0045] This step utilizes the parallel particle swarm optimization algorithm, a specific choice tailored to address the pain points of coalfield hydrological model calibration, rather than a simple application of general algorithms. Its core advantages are as follows: Global optimization capability, adaptable to high-dimensional nonlinear problems: Particle swarm optimization (PSO) is a typical swarm intelligence global optimization algorithm. Simulating bird flock foraging behavior, it uses multiple particles to search in parallel within the parameter space, with each particle corresponding to a complete set of hydrological parameters to be optimized. By iteratively updating the individual optimal and global optimal positions, it achieves global optimization in high-dimensional parameter space, perfectly adapting to the synchronous optimization of dozens or even hundreds of hydrological parameters in coalfield models, fundamentally avoiding the local optima problem of manual trial and error. Parallel computing architecture, overcoming efficiency bottlenecks: Addressing the pain points of large computational load and numerous iterations in a single simulation of large-scale 3D coalfield models, this algorithm employs a parallel computing architecture. The simulation computation tasks of dozens of particles are simultaneously distributed to multiple CPU cores / computing nodes for parallel execution. This compresses the calibration work of traditional single-threaded algorithms, which takes weeks, to complete in hours to days, completely solving the efficiency bottleneck of multi-parameter global optimization and enabling engineering implementation. Quantifying the objective function to ensure global fitting accuracy: This step clearly defines the minimum root mean square error (RMSE) between the simulated water level and the measured water level as the sole optimization objective. RMSE is a recognized absolute deviation quantification index in the field of hydrological models. It is extremely sensitive to large deviation values and can force the model to avoid the problem of large deviations in water levels in local areas, ensuring the fitting accuracy of the model in the entire study area and throughout the entire time period, rather than just fitting the water level values of a few observation wells.
[0046] Parameter space delineation under geological constraints: For each hydrogeological parameter to be calibrated, based on coalfield pumping tests, regional hydrogeological manuals, and measured data from adjacent areas, set the upper and lower limits of the parameters that conform to geological laws, and strictly prohibit extreme values without physical meaning in parameter optimization, so as to avoid invalid optimization of randomly adjusting parameters in order to fit the water level from the root cause. Particle swarm initialization: Randomly generate 30 - 100 particles within the delineated parameter space. Each particle corresponds to a complete set of hydrogeological parameter combinations to form an initial search population. Parallel simulation and fitness calculation: Input the parameter combinations of each particle into the initial numerical model, and perform parallel numerical simulations of groundwater flow during the calibration period to calculate the RMSE value between the simulated water level and the measured water level corresponding to each particle as the fitness value of the particle (the smaller the RMSE, the higher the fitness). Iterative update and optimization: Update the individual optimal position of each particle and the global optimal position of the entire population according to the fitness value, adjust the flight speed and direction of the particles, generate new parameter combinations, and repeat parallel simulation and fitness calculation. Convergence termination: When the number of iterations reaches the preset upper limit, or the global optimal RMSE value has not decreased significantly for more than 20 generations, stop the iteration, output the globally optimal hydrogeological parameter combination, and complete the automatic calibration.
[0047] This step is the final access threshold for the benchmark numerical simulation model, establishing a strict quantitative qualification standard for the coalfield hydrogeological model, and completely solving the problems of no standard and strong subjectivity in traditional model verification.
[0048] Select the Nash efficiency coefficient (NSE) as the core verification index, which is a recognized gold standard in the field of hydrogeological models, and forms a perfect complement with the RMSE used during the calibration period: RMSE is an absolute deviation index used for global fitting optimization during the calibration period; NSE is a relative fitting accuracy index specifically used to evaluate the prediction ability of the model, and its calculation formula is: .
[0049] Among them, , is the measured water level, , is the simulated water level, , is the average value of the measured water level, and n is the number of observation data during the verification period; The value range of NSE is (−∞, 1]. The closer it is to 1, the stronger the prediction ability of the model: when NSE = 1, the simulated value is exactly the same as the measured value; when 0.75 < NSE < 1, the model accuracy is excellent; when 0.65 < NSE ≤ 0.75, the model accuracy is good; when NSE ≤ 0.5, the model results are completely untrustworthy and have no predictive value.
[0050] In response to the high safety requirements of coal mine water control, this invention sets the preset threshold as NSE≥0.75 (excellent level), which is much higher than the industry-standard 0.65. This is because the results of coal mine groundwater simulation are directly related to the safe production of underground water inrush prevention and control, and the requirements for prediction accuracy are much higher than those of ordinary hydrogeological surveys.
[0051] The core execution logic of the validation is as follows: Keeping the optimal hydrological parameters completely fixed after calibration, the source-sink terms and boundary conditions for the validation period are input into the model. Groundwater flow simulations for the validation period are run, and simulated water level time-series data from each observation well are output. The NSE (Numerical Sequence) values between the simulated and measured values during the validation period are calculated, and the RMSE (Real-Time Sequence) and coefficient of determination R² (R²) are simultaneously validated to form a multi-dimensional accuracy evaluation system. If the NSE reaches the preset threshold, it indicates that the model can not only fit the known data from the calibration period but also accurately predict groundwater dynamics under independent time periods and conditions, possessing stable generalization prediction capabilities, and is formally determined as the benchmark numerical simulation model. If the NSE does not reach the threshold, it indicates that the model has overfitting issues or unreasonable geological structure and boundary condition settings. The initial model needs to be optimized, and calibration and validation re-executed until the accuracy requirements are met.
[0052] In this embodiment, the dynamic hydrological parameter field is obtained in step S3 by: denoising and normalizing the temperature and strain time series data to extract their spatiotemporal variation characteristics; constructing an inversion model based on the heat conduction equation and the thermo-solid coupling constitutive relation, using the preprocessed temperature and strain time series data as input to the inversion model, and using the aquifer permeability coefficient and porosity as the parameters to be inverted; using a Bayesian inversion framework and combining the Markov chain Monte Carlo method to solve the inversion model to obtain the posterior probability distribution of the parameters to be inverted; and extracting the mean of the posterior probability distribution as the dynamic hydrological parameter field.
[0053] It should be noted that the temperature and strain time-series data processed comes from a specially deployed distributed fiber optic sensing system (DTS distributed temperature sensing + DAS distributed strain sensing) in the target coalfield aquifer. This is fundamentally different from traditional point sensors and borehole logging data: it is spatially continuous and temporally high-frequency monitoring data for the entire well section, with a spatial resolution of 0.5-1m and a temporal sampling frequency of minutes. It can completely capture the dynamic response of temperature and strain in the entire vertical section and throughout the entire time period of the aquifer. It is currently the only technical means in the coalfield field that can achieve continuous real-time perception of the entire aquifer seepage activity.
[0054] Coalfield underground fiber optic monitoring faces extremely complex conditions such as strong vibrations during mining, underground electromagnetic interference, fiber creep, geothermal field background drift, and heat dissipation interference from mining. The preprocessing in this step is specifically designed for this scenario, rather than general data processing: Targeted denoising: For DTS temperature data, wavelet adaptive threshold denoising is used to remove peak noise caused by mining vibration and electromagnetic interference. At the same time, sliding window regression is used to eliminate static background drift caused by geothermal gradient and heat dissipation from underground equipment. For DAS strain data, singular value decomposition + moving average filtering is used to eliminate low-frequency trend drift caused by fiber installation stress release and rock creep. At the same time, pulse-like abnormal noise caused by blasting at the working face and equipment start-up and shutdown is removed. Spatiotemporal registration and normalization: First, complete the spatiotemporal alignment of DTS and DAS data, unify the spatial sampling interval and time step of the two types of data, and eliminate the sampling misalignment of different sensing units; then, through min-max normalization, map the two physical quantities of temperature and strain with different dimensions to the [0,1] interval, eliminate the weight interference of the difference in dimensions on the subsequent inversion, and ensure that the constraint effect of the two physical fields is balanced and effective.
[0055] The extracted features are not statistical characteristics of the original data, but dynamic response features directly related to groundwater seepage activity. These core features include: the temperature change rate and strain dynamic response rate in the time dimension, and the temperature anomaly gradient and strain abrupt change segments in the spatial dimension. By eliminating static background fields and anthropogenic interference features unrelated to seepage, only effective dynamic features caused by groundwater flow and pore water pressure changes in the aquifer are retained. This fundamentally reduces invalid inputs to the inversion model and lowers the non-uniqueness of the inversion.
[0056] Traditional distributed fiber optic hydrological monitoring in coalfields relies solely on DTS single-temperature field data and uses pure heat conduction-convection equations to invert permeability coefficients. This approach faces two insurmountable industry bottlenecks: First, the inversion is highly non-unique: temperature changes can be caused by groundwater seepage and convection, as well as by variations in geothermal gradients, mining heat dissipation, and fluctuations in underground environmental temperature. A single field constraint cannot distinguish between interference signals and valid seepage signals. Second, the parameter dimensions are severely insufficient: only permeability coefficients can be inverted, failing to simultaneously acquire porosity. Porosity is a core parameter determining aquifer storage capacity, rock compressibility, and pore water pressure propagation patterns, and is an indispensable key input for groundwater numerical simulation. Single-parameter inversion cannot provide complete parameter constraints for the numerical model.
[0057] The constructed inversion model is a dedicated inversion framework that fully couples the temperature field, stress field, and seepage field, rather than a simple superposition of heat conduction and solid mechanics. Its core is to lock in a unique mapping relationship between the parameters to be inverted and the monitoring data through the dual constraints of two physical fields. The core logic is as follows: Temperature field governing equation (heat conduction-convective coupling): Based on the convective heat transfer of groundwater seepage, a quantitative relationship between temperature changes and seepage parameters is established. The strength of the convective heat transfer term in the equation is directly determined by the aquifer permeability coefficient (which determines the Darcy flow velocity) and porosity (which determines the groundwater volume fraction), establishing a direct correlation between temperature data and the two parameters to be inverted. Thermal-solid coupling constitutive relationship (stress-seepage coupling): Based on Terzaghi's effective stress principle, a quantitative relationship between strain changes and seepage parameters is established. The strain response of the aquifer rock mass originates from two aspects: first, the strain caused by changes in pore water pressure due to groundwater seepage, which in turn leads to changes in the effective stress of the rock mass; second, the thermal expansion and contraction strain of the rock mass caused by temperature changes. The propagation velocity of pore water pressure is determined by the permeability coefficient, and the compressive deformation characteristics of the rock mass are determined by porosity. A second independent constraint is established between the strain data and the two parameters to be inverted. The forward-inversion closed-loop logic is as follows: the forward operator of the inversion model is the three-field fully coupled control equation mentioned above. By inputting a set of parameters to be inverted (permeability coefficient, porosity), the corresponding spatiotemporal temperature and strain time-series data can be calculated using forward modeling. The core objective of the inversion is to minimize the residual between the forward calculation results and the fiber optic measured data.
[0058] This inversion model has been customized and optimized specifically for coalfield mining conditions: dynamic boundary conditions such as working face advancement, mining and unloading disturbance, and mine drainage intensity are incorporated into the model, rather than using static boundaries under ideal laboratory conditions. This ensures that the inversion model can accurately capture the dynamic changes of aquifer parameters under mining disturbance and adapt to the real-world scenario of safe coal mine production.
[0059] Traditional deterministic inversion methods (such as least squares) can only output a single set of optimal parameter values, which has two fatal flaws: first, it is prone to getting trapped in local optima and cannot find the globally optimal parameter distribution; second, it cannot quantify the uncertainty of the parameters, cannot determine whether the output optimal value has engineering reliability, and is very likely to cause systematic deviations in subsequent numerical simulations. This section fundamentally solves the above problems through a Bayesian probabilistic inversion framework.
[0060] Based on Bayes' theorem, this step defines the permeability coefficient and porosity to be inverted as random variables. Combining prior information and measured data, the posterior probability distribution of the parameters is solved. The core formula is: Where: θ represents the parameters to be inverted (permeability coefficient, porosity), and D represents the preprocessed measured data of fiber temperature and strain. The prior probability distribution is not an unconstrained random distribution, but a parameter distribution that conforms to geological laws, based on the hydrological parameters of the benchmark model that has been calibrated and the hydrogeological test results of the coalfield area (such as the permeability coefficient adopting the log-normal distribution recognized in the field of hydrogeology). This not only reduces the parameter search space and improves the inversion convergence efficiency, but also ensures the physical compatibility between the inversion parameters and the benchmark model. The likelihood function is used to measure the degree of matching between the forward modeling results and the measured data. It is constructed using a multivariate normal distribution and incorporates the residuals of the two physical fields, temperature and strain. Appropriate weights are assigned to the two fields to ensure that the dual constraints are balanced and effective, and to avoid a single field dominating the inversion process. The posterior probability distribution to be determined represents all possible values of the parameter to be inverted and their corresponding probabilities after combining all measured data and prior information. It is the core objective of the inversion solution.
[0061] This step employs the adaptive DRAM-MCMC algorithm to solve for the posterior probability distribution, a specific choice for large-scale 3D inversion scenarios in coalfields. The core algorithm logic involves generating a Markov chain in a high-dimensional parameter space for random sampling. When the Markov chain converges, the collected samples conform to the posterior probability distribution of the parameters to be inverted. The adaptive optimization design addresses the shortcomings of traditional MCMC algorithms, such as slow convergence and low sampling efficiency in high-dimensional spaces, by adaptively adjusting the sampling step size and proposal distribution. This efficiently solves for the posterior distribution of parameters across the entire 3D space of the coalfield. Furthermore, it not only obtains the optimal parameter values but also quantifies the uncertainty of parameters at each spatial location, providing confident parameter constraints for subsequent PINN parameter optimization. This avoids unreliable parameter values misleading the optimization process and ensures the physical reliability of the entire parameter optimization process.
[0062] The core reason for choosing the mean of the posterior probability distribution as the final parameter value, rather than the single-point maximum a posteriori (MAP) value, is that the posterior mean is the mathematical expectation of all convergent samples, which integrates all possible values of the parameter. It is more robust than the single-point MAP value, can eliminate the bias caused by sampling randomness and inversion uncertainty to the greatest extent, and has stronger compatibility with the parameter field of the benchmark model, making it more suitable for the iterative needs of subsequent numerical simulation and parameter optimization.
[0063] The output dynamic hydrological parameter field possesses three core characteristics that traditional methods cannot achieve: High spatial resolution: The spatial resolution of the parameter field is consistent with that of fiber optic sensing (0.5-1m), realizing continuous parameter distribution along the entire vertical section of the aquifer and along the fiber optic deployment path in the plane. This completely fills the parameter gap of hundreds of meters between traditional boreholes, providing unprecedentedly fine parameter constraints for three-dimensional numerical simulation; Dynamic time update: The parameter field is not a one-time static result, but a time-series parameter sequence that is updated in real time with fiber optic monitoring. Every fixed time step (e.g., 1-4 hours), the system collects a new set of temperature and strain data, repeats the above preprocessing-inversion-solution process, and can generate the hydrological parameter field at the corresponding time, accurately capturing the dynamic changes in permeability and porosity of the aquifer caused by mining unloading and fracture development during coal mining; Full three-dimensional spatial coverage: Through the network deployment of distributed optical fibers in multiple boreholes, combined with spatial interpolation algorithms, a dynamic three-dimensional hydrological parameter field of the entire target coalfield can be generated, rather than a one-dimensional parameter sequence in the vertical direction of a single well, which is fully adaptable to the needs of subsequent three-dimensional groundwater numerical simulation.
[0064] In this embodiment, step S4 involves solving the loss function embedded in the partial differential equation of groundwater flow. Specifically, this involves: constructing a deep neural network consisting of an input layer, multiple fully connected hidden layers, and an output layer. The input layer receives spatial and temporal coordinates, and the output layer outputs the predicted hydraulic head value. A total loss function is constructed, including data fitting terms, initial condition constraints, boundary condition constraints, and physical equation constraints. The physical equation constraints are based on the groundwater flow continuity equation and Darcy's law, and are used to constrain the predicted hydraulic head value output by the deep neural network to satisfy the physical laws of groundwater flow. An adaptive moment estimation optimizer is used to iteratively update the network weights of the deep neural network with the goal of minimizing the total loss function. After training, optimized hydrological parameters are extracted from the deep neural network.
[0065] It is important to note that the neural network constructed in this section is essentially a global approximator of the spatiotemporal continuous function of groundwater head. This is a fundamental difference from traditional neural network applications in the hydrological field: traditional hydrological neural networks primarily take discrete features such as monitoring data and borehole parameters as inputs, and output hydrological parameters or single-point water levels, only capable of fitting discrete points and unable to characterize the continuous seepage patterns across the entire aquifer. This architecture uses four-dimensional spatiotemporal coordinates (x, y, z, t) as the sole input, where (x, y, z) are the three-dimensional spatial coordinates of the aquifer, t is the time coordinate, and the output is the continuously predicted head value h(x, y, z, t) at the corresponding spatiotemporal location. Its core principle is that groundwater head is a continuous function satisfying partial differential equation constraints in the spatiotemporal domain, and fully connected deep neural networks possess universal approximation capabilities, accurately fitting any nonlinear continuous function. This fundamentally solves the accuracy bottleneck caused by the discrete grids of traditional numerical simulations and the deficiency of traditional neural networks in achieving continuous prediction across the entire domain.
[0066] This architecture is not a general fully connected network, but is specifically optimized for the geological and engineering characteristics of coalfield aquifers: Hidden layer adaptation design: To address the strong heterogeneity of vertical stratification and planar partitioning in coalfield aquifers, 6-8 fully connected hidden layers are used, with 50-100 neurons in each layer, balancing nonlinear fitting capability and training convergence efficiency; The activation function is specifically chosen to be the Tanh or Swish function, rather than the general ReLU function. The core reason is that PINN needs to solve the higher-order partial derivatives of the hydraulic head with respect to the spatiotemporal coordinates through automatic differentiation, while the higher-order derivatives of the ReLU function are discontinuous, which will lead to distortion of the physical equation residual calculation. The Tanh / Swish function has the characteristic of infinite differentiability, which is fully adapted to the solution requirements of the partial differential equation of groundwater flow. Spatiotemporal adaptability: Considering the long temporal and large spatial scale characteristics of coalfield mining, the spatiotemporal coordinates of the network input are normalized and decoupled, mapping spatial and temporal coordinates to the [0,1] interval respectively. This avoids the gradient vanishing problem caused by large spatial scales and long temporal periods, ensuring prediction accuracy throughout the entire cycle from the initial mining phase to stable drainage. Parameter coupling design: This network sets the core hydrological parameters to be optimized (permeability coefficient tensor, porosity, and water storage coefficient) as learnable parameters trained synchronously with the network weights and biases, rather than fixed input or output values. This is the most critical proprietary design of this architecture: it fully couples hydrological parameter optimization with head calculation, simultaneously completing the solution of the global continuous head field and the global optimization of hydrological parameters during network training.
[0067] This step is the core distinguishing feature of PINN from purely data-driven neural networks, and it is also the core carrier for this invention to achieve dual data-physical driving. It constructs a dual closed loop of measured data constraints and hydrogeological law constraints through a four-component weighted loss function, fundamentally solving the fatal flaws of purely data-driven models, such as black-box nature, lack of physical meaning, overfitting, and poor generalization.
[0068] The total loss function is a weighted sum of the four components, and the formula is: Where ω1−ω4 are the weight coefficients of each loss term, and the specific design logic for the four loss terms in the coalfield scenario is as follows: Data fitting term The core principle is to constrain the neural network to minimize the deviation between the predicted water head and the measured true value at spatiotemporal points with available measured data, using mean squared error (MSE) for construction. Its input measured data includes three distinct core data types: measured water level time-series data from observation wells, simulated water level data from benchmark numerical simulation models, and water head distribution data corresponding to the dynamic hydrological parameter field retrieved via fiber optic inversion; initial condition constraints... The output head of the constrained neural network at the initial time (t=0) perfectly matches the initial head field after calibration of the benchmark model, and is also constructed using MSE. Coal mine groundwater simulation is a long-term dynamic process; the impact of face mining and drainage accumulates over time, and even small deviations in initial conditions can lead to an exponential amplification of errors in subsequent long-term predictions. This invention locks in the physical reality of the initial flow field from the training start point, completely solving the industry problem of large initial time deviations and cumulative divergence of long-term prediction errors in traditional time-series neural networks. The core challenge of coalfield hydrological monitoring is the extremely limited number and highly uneven spatial distribution of observation wells. Traditional pure data-driven models, relying solely on this factor for fitting, can only guarantee the fitting effect around the observation wells; the entire region without observation data is completely distorted. In this invention, this factor is only one constraint term in the total loss function, not the sole objective, ensuring both the fitting accuracy of measured data and preventing the failure of global predictions due to data sparsity; boundary condition constraint terms... The constrained neural network outputs hydraulic head values at spatiotemporal points along aquifer boundaries that strictly satisfy the specific hydrogeological boundary conditions of coalfields, including constant hydraulic head boundaries (river recharge, aquifer outcrops), constant flow boundaries (mine drainage, working face dewatering), and impermeable boundaries (closed faults, stratigraphic pinch-outs), all constructed using the MSE (Mean Segregated Segregated Equations). The boundary conditions of coalfield aquifers determine the core patterns of groundwater recharge, runoff, and discharge throughout the region, and are crucial for flow field simulation. Traditional pure data-driven models cannot handle hydrogeological boundaries, leading to completely distorted predictions in the boundary region and even results that violate basic hydrological laws (such as sudden flow changes at impermeable boundaries). This forced constraint network output conforms to the physical rules of the boundaries, ensuring the overall rationality of the flow field across the entire region; the physical equation constraint terms... This is the core innovation of PINN and a landmark design that distinguishes this invention from all traditional coalfield hydrological parameter optimization methods. Its core is to embed the three-dimensional unsteady groundwater flow control equations (based on the continuity equation and Darcy's law) into the loss function. Through the automatic differentiation mechanism of the neural network, the residuals of the control equations are calculated, thereby constraining the network output to strictly follow the fundamental physical laws of groundwater flow.
[0069] The governing equations for the three-dimensional unsteady flow in coalfield aquifers are: in, The water storage coefficient, Let h be the permeability coefficient in three directions, h be the predicted head, and W be the source and sink terms (mine drainage, recharge, etc.). The left side of the equation represents the head variation over time, and the right side represents the seepage and source / sink terms. It is entirely based on Darcy's law and the law of conservation of mass, representing the core physical laws of groundwater flow. Using the automatic differentiation mechanism of a neural network, the first and second partial derivatives of the predicted head h output by the network are calculated with respect to the input spatial coordinates (x, y, z), and the first partial derivative is calculated with respect to the time coordinate t. The obtained partial derivatives and the hydrological parameters to be optimized are substituted into the above governing equations, and the residuals (i.e., the difference between the left and right sides of the equation) are calculated. The MSE value of the residuals is... .
[0070] This is equivalent to setting a virtual physical constraint point at every spatiotemporal coordinate point across the entire aquifer. Even in areas without any actual monitoring data, the neural network output must strictly satisfy the physical laws of groundwater flow. This completely solves the core industry dilemma of sparse coalfield observation wells and the inability to constrain data-free areas, upgrading parameter optimization from merely fitting a few observation points to conforming to the hydrogeological and physical laws across the entire region. It fundamentally avoids the industry's persistent problem of adjusting parameters to have no physical meaning in order to fit water levels.
[0071] In this embodiment, step S5 involves multi-scale simulation calculations, specifically: constructing a nested grid system, using locally refined grids in the well periphery and geologically complex areas, and coarse grids in other areas; mapping the optimized hydrological parameter set onto each grid cell of the nested grid system; discretizing the groundwater flow control equations using the finite volume method based on the nested grid system to construct a large-scale sparse linear equation set; and using the algebraic multigrid preprocessing conjugate gradient method to solve the large-scale sparse linear equation set in parallel, obtaining high-precision prediction results for aquifer water level and flow field distribution.
[0072] It should be noted that the constructed nested grid system is designed entirely around the core engineering needs of overall control of the flow field in the entire coalfield for water prevention and control, and precise focus on key risk areas, completely breaking the deadlock that traditional uniform grids cannot achieve both accuracy and efficiency.
[0073] Traditional numerical simulations of coalfields generally use uniform grids across the entire area. If a coarse grid (such as 50m×50m) is used, it will completely smooth out the drastic changes in the flow field in key areas such as the well perimeter and fault zones, making it impossible to capture potential water-conducting channels and water inrush risk points. The prediction results will not have any guiding value for water prevention and control. If a fine grid (such as 5m×5m) is used across the entire area, the number of grids in a three-dimensional model of a coalfield covering tens of square kilometers will exceed hundreds of millions, and a single simulation calculation will take several days, which is completely unsuitable for the needs of dynamic real-time prediction of coal mining.
[0074] This nested grid system adopts a two-level seamless nested architecture of a global background coarse grid and multi-segment locally refined sub-grids, precisely adapted to the risk characteristics of coalfields: Global background coarse grid: Covers the entire target coalfield area, with a grid size typically set to 50m×50m, vertically matching the layer thickness of aquifers / impermeable layers. Its core function is to control the overall recharge-runoff-discharge pattern of groundwater throughout the coalfield, ensuring the computational integrity of the global flow field while keeping the basic computational workload within a reasonable range. Locally refined grid: Grid refinement is applied only to high-risk areas of core concern for water control in two types of coal mines, with refinement sizes reaching 5m×5m or even 1m×1m. It achieves seamless boundary nesting with the global coarse grid and complete flow conservation: Areas surrounding wells: Includes hydrological observation wells, mine drainage wells, longwall face roadways, and areas surrounding tunneling heads. These areas exhibit the most dramatic changes in groundwater head gradient and flow velocity, and are also core areas for mine drainage and water hazard monitoring. Densified grids can accurately capture the water level drawdown cones and abrupt changes in seepage velocity around the mine shafts, avoiding the head calculation errors caused by coarse grids. Complex geological structures include faults, collapse columns, fracture zones, lithological facies transition zones, and aquifer outcrops. These areas are the core sources of underground water inrush in coal mines. Traditional coarse grids completely obscure the strong water-conducting characteristics of structural zones, failing to identify potential water-conducting channels. Densified grids can accurately depict the spatial morphology and parameter abrupt changes of structural zones, restoring their true seepage patterns. Vertically tailored: For the geological characteristics of layered aquifers in coalfields, vertical grid densification is applied to key water-bearing aquifers on the roof and floor of the main coal seam, while coarser grids are used for weakly permeable aquitards. This ensures the vertical seepage accuracy of key aquifers while avoiding unnecessary computational waste.
[0075] It solves the fundamental contradiction of insufficient accuracy of traditional uniform grid coarse grids and explosive computational load of fine grids. Under the premise of controlling the overall computational load, it achieves ultra-high precision characterization of the core risk area of water control, so that the numerical simulation results can truly have practical value in guiding on-site water control projects.
[0076] The entire continuous hydrological parameter field output in the previous step is completely and losslessly transformed into parameter values for each discrete grid cell that can be calculated in numerical simulation. This completely solves the industry problem of losing high-resolution parameter information during discretization and the inability to implement the advantages of optimized parameters in the traditional parameter assignment process.
[0077] This mapping is not a simple spatial interpolation, but a dedicated adaptation scheme for nested grids and continuous parameter fields. The core design is as follows: Integral-average mapping principle: The hydrological parameter set obtained from prior optimization is a continuous function corresponding to three-dimensional spatial coordinates (not discrete point parameter values). For each cell in the nested grid, the integral average value of the continuous parameter field within the spatial range of that cell is first calculated during mapping, and then assigned to that grid cell, instead of using single-point interpolation at the cell center point. This design ensures that regardless of whether the grid is refined or coarsened, the discretized parameter field can completely restore the overall seepage characteristics of the original continuous parameter field, without parameter distortion due to changes in grid scale. Complete preservation of anisotropic parameters: Given the strong anisotropy of coalfield layered rock strata, where the difference in horizontal and vertical permeability coefficients can reach 2-4 orders of magnitude, the mapping process completely assigns the optimized permeability coefficient tensor (permeability coefficients in the x, y, and z directions) to each grid cell, instead of using traditional equivalent homogeneous parameters, accurately restoring the bedding-parameter and vertical seepage differences of layered aquifers. Parameter assignment adapted to grid scale: For fine grid cells, the high-resolution heterogeneity of the continuous parameter field is strictly matched to accurately capture the parameter abrupt change characteristics of fault zones and lithological facies transition zones; for coarse grid cells across the entire domain, equivalent seepage parameters are used for calculation to ensure that the equivalent permeability of the coarse grid is completely consistent with the overall seepage characteristics of the fine parameter field in the corresponding region, avoiding distortion of seepage patterns caused by coarsening.
[0078] Traditional coalfield models typically employ layered homogeneous parameter assignment or discrete borehole point interpolation, resulting in the complete loss of high-resolution, highly heterogeneous parameter information obtained from previous optimization during the discretization process. Consequently, the technical advantages of parameter optimization cannot be translated into improved simulation accuracy. This step utilizes a dedicated mapping scheme to fully transfer the optimized parameter accuracy to the simulation mesh, fundamentally ensuring the physical realism of the numerical simulation.
[0079] The finite volume method (FVM) was specifically chosen for equation discretization. The core reason for this is to adapt to the rigid engineering requirements of strict water conservation in coalfield water control. This is the most fundamental difference between FVM and traditional finite difference method (FDM) and finite element method (FEM).
[0080] In coal mine water control, the prediction of mine drainage volume, calculation of inrush flow, and assessment of aquifer dewatering volume must strictly adhere to water balance (mass conservation). Otherwise, the predicted water level and flow rate results will exhibit systematic deviations, rendering them completely unsuitable for guiding on-site engineering. The core advantage of the finite volume method lies in the dual mass conservation of both local cells and the global computational domain: it performs volume integration on the groundwater flow control equations in each grid cell, ensuring that the inflow rate = outflow rate + change in storage within the cell for each cell, thus eliminating the problem of water imbalance from the perspective of discretization principles.
[0081] For the nested grid system of this invention, the finite volume method was specifically used for discretization optimization, which fundamentally solves the conservation problem of nested grid boundaries: A flux-matched discretization scheme is adopted for the nested boundaries of coarse and fine sub-grids, strictly ensuring the continuity of flow and head between the coarse and fine grids at the boundaries. There is neither spurious water replenishment nor water loss, completely solving the industry-wide problem of boundary non-conservation and flow field distortion during traditional nested grid discretization. The discretized governing equations are completely consistent with the three-dimensional unsteady groundwater flow governing equations used in the previous baseline model and parameter optimization stage, ensuring the self-consistency of the physical laws throughout the entire invention process and avoiding systematic deviations caused by equation differences. After discretization, a large-scale sparse linear equation system is finally formed: each grid cell is hydraulically connected only to its adjacent cells, and the vast majority of elements in the coefficient matrix of the equation system are 0, exhibiting typical sparsity characteristics, providing a foundation for subsequent efficient parallel solutions.
[0082] Traditional finite difference methods cannot guarantee local water conservation in non-uniform, nested mesh environments. The finite element method suffers from high computational cost and poor water conservation, leading to water balance errors exceeding 30% in large-scale coalfield models. Consequently, the predicted mine drainage and water level changes often deviate significantly from actual field conditions. This approach utilizes the conservation discretization of the finite volume method to fundamentally ensure water balance in the simulation results, giving the predictions practical engineering value.
[0083] The conjugate gradient method (CG) with algebraic multigrid (AMG) preprocessing is specifically selected, which completely solves the industry bottleneck of slow convergence, extremely low computational efficiency, and inability to meet the real-time prediction requirements of large-scale nested grid models in coalfields.
[0084] The linear equations discretized from the nested grid model of the coalfield can have hundreds of thousands or even tens of millions of unknowns. Furthermore, due to the scale difference of more than 50 times between the coarse and fine grids, the condition number of the coefficient matrix is extremely poor. Traditional solvers such as Gaussian elimination and Jacobi iteration are either completely unable to solve the problem or have extremely slow convergence speeds. A single simulation can take up to several days, which is completely unsuitable for the real-time prediction needs of daily progress and dynamic changes in the flow field of the coal mine working face.
[0085] This solution scheme consists of two core modules, deeply adapted to the characteristics of coalfield models: The first is the Conjugate Gradient (CG) core solver. The coefficient matrix of the groundwater flow control equations, after discretization, is a typical symmetric positive definite matrix, perfectly suited to the requirements of the CG method. Compared to traditional iterative methods, the CG method improves convergence speed by more than an order of magnitude with extremely low memory usage, perfectly adapting to solving large-scale sparse equation systems. The second is the Algebraic Multigrid (AMG) preprocessing module. This is the core module for solving the convergence problem of nested grids. Addressing the issue of extremely poor condition numbers in the coefficient matrix and convergence stagnation in the CG method due to large differences in grid scales, the AMG preprocessing module constructs multi-layer grid constraints and interpolation operators to smooth the residuals of the equation system across grids of different scales. This reduces the condition number of the coefficient matrix by 2-3 orders of magnitude, further improving the convergence speed of the CG method by more than an order of magnitude, completely solving the convergence problem of nested grid models. Parallel optimization by domain decomposition: For the large-scale computational domain of coalfields, a domain decomposition strategy is adopted to divide the entire grid system into multiple sub-regions. Each sub-region is assigned to an independent CPU core / computing node for parallel solution. The boundary information between sub-regions is synchronized through MPI communication, which can achieve a near-linear parallel speedup ratio. The single simulation calculation that originally required several days is compressed to minutes, which fully meets the engineering requirements of real-time dynamic prediction in coal mines.
[0086] After the solution is completed, not only is the head value of each grid cell obtained, but also the full flow field parameters such as seepage velocity, hydraulic gradient, and unit width flow rate of each cell are calculated simultaneously based on Darcy's law. Finally, a complete prediction result of aquifer water level and flow field distribution with high spatiotemporal resolution is formed, which directly provides core data support for the subsequent identification of water-rich anomaly areas and potential water-conducting channels.
[0087] In this embodiment, step S6 is also included: performing three-dimensional visualization rendering of the high-precision aquifer water level and flow field distribution prediction results to generate a dynamic interactive display interface; based on the dynamic interactive display interface, identifying and delineating water-rich anomaly areas and potential water-conducting channels; and outputting the boundary information of water-rich anomaly areas and potential water-conducting channels as the target area coordinate set for coal mine water control work.
[0088] It should be noted that the output of traditional numerical simulations of groundwater in coalfields is generally in the form of two-dimensional water level contour maps, static three-dimensional slice maps, or discrete grid data within specialized software. These results can only be interpreted by hydrogeological professionals. Production managers and underground construction workers at the coal mine site cannot quickly access the core risk information, resulting in a large amount of high-precision simulation results being locked in specialized computers and unable to be translated into practical water control measures on-site. The core objective of this section is to break down these professional barriers and transform the obscure discrete numerical results into a visually interactive tool that is understandable, operable, and decision-making-friendly for all positions.
[0089] The rendering system does not display flow field data in isolation, but deeply integrates all core elements of coalfield water control with flow field results to construct a 1:1 3D scene that restores the real geological and engineering conditions of the mining area: Geological elements: It fully integrates the spatial morphology of core geological bodies such as aquifers, aquitards, main coal seams, faults, collapse columns, and lithological facies transition zones, achieving precise spatial registration with flow field data; Engineering elements: It simultaneously embeds underground engineering facilities such as longwall faces, tunnels, drainage wells, observation wells, and advanced drilling holes, intuitively showing the spatial relationship between water hazard risk sources and mining engineering; Flow field-specific rendering: It uses volume rendering technology to display the continuous 3D spatial distribution of aquifer water-bearing capacity, isosurface rendering to display iso-water levels and hydraulic gradient critical surfaces, and the color mapping system is fully adapted to the cognitive habits of coal mine sites (high-risk water-bearing areas and strong seepage areas are highlighted with red tones, and low-risk areas are highlighted with blue tones), allowing for rapid identification of risk distribution without professional knowledge.
[0090] The interactive interface is developed entirely around the on-site decision-making needs of coal mine water control. Core functions are all industry-specific designs, rather than generic zoom and rotate operations: Spatiotemporal dynamic interaction: By sliding along the timeline, the dynamic evolution of aquifer water level, flow field, and water-bearing capacity can be displayed in real time throughout the entire cycle of face advancement and mine drainage. It can also predict water hazard risk changes under the next 7-day and 30-day mining progress, completely solving the pain point of traditional static results being unable to predict dynamic risks; Engineering-based cross-section interaction: Supports arbitrary planar cross-sections, vertical cross-sections, and directional cross-sections along the roadway, allowing direct viewing of the area in front of the tunnel face and the mining area. The system accurately obtains vertical risk information directly relevant to underground construction by analyzing the aquifer's water content and hydraulic gradient distribution on the top and bottom plates of the working face; it features one-click parameter query interaction: clicking on any spatial location within the 3D scene will bring up core parameters such as the mining area coordinates, aquifer thickness, permeability coefficient, water level, hydraulic gradient, and seepage velocity, allowing even non-professionals to quickly obtain accurate quantitative data; and it supports customized early warning interaction: allowing on-site personnel to customize safety thresholds for water level, hydraulic gradient, and seepage velocity according to mining area water prevention and control regulations. The interface automatically highlights high-risk areas exceeding the thresholds, directly connecting to the on-site water hazard early warning system.
[0091] It breaks down the professional barriers of coalfield hydrological simulation results, transforming numerical results that could only be interpreted by hydrogeological experts into decision-making tools that can be understood and used by everyone from the management to the underground construction teams in coal mines, thus fundamentally solving the industry problem of the disconnect between simulation and on-site operations.
[0092] A human-machine collaborative identification system combining algorithm pre-identification and human experience verification was constructed, perfectly balancing the objectivity of the algorithm with the experience of on-site engineers. This system solves two major defects of traditional identification methods: First, purely manual interpretation of contour maps is highly subjective, inefficient, and inaccurate in judging three-dimensional spatial relationships, making it easy to miss deep water-conducting channels; second, purely automatic algorithm identification does not take into account the actual geological conditions of the mining area, mining history, and water hazard cases, resulting in a high rate of misjudgments and false anomalies, and the delineated target area is seriously inconsistent with the actual situation on site.
[0093] The identification and delineation in this stage are completed entirely within a dynamic interactive interface, forming a closed-loop workflow. The core steps are: Algorithm pre-identification and highlighting: Based on high-precision flow field simulation results, the system automatically completes the initial screening of risks across the entire area, highlighting potential risk areas that meet the characteristics of water-rich anomalies and water-conducting channels in the 3D interactive interface, providing targeted guidance for manual interpretation and significantly reducing the workload of manual investigation; Multi-dimensional manual verification: Hydrogeologists and water control personnel at the coal mine site, combined with their on-site experience in the geological patterns of the mining area, fault water conductivity, borehole exposure, and historical water hazard cases, verify the information through interface segmentation and timelines. Dynamic simulation and parameter query functions are used to verify the pre-identified areas one by one: false anomalies caused by local parameter fluctuations are eliminated, hidden risk areas that were not identified by the algorithm but conform to the geological laws on site are supplemented, and the boundary range of the risk areas is corrected; precise three-dimensional spatial delineation: the final delineation is not a two-dimensional plane range, but a closed and continuous geological body boundary in three-dimensional space, clearly marking the top and bottom plate depths and planar distribution range of the water-rich anomaly area, the spatial direction and extension scale of the potential water diversion channel, as well as its core engineering parameters such as vertical and horizontal distances with the mining face and roadway, which fully match the design requirements of underground water control engineering.
[0094] It solves the industry dilemma of traditional water hazard risk identification being either purely subjective or purely mechanical. It avoids omissions and misjudgments caused by manual interpretation and makes up for the shortcomings of pure algorithm identification that are detached from the actual situation on site. It significantly improves the accuracy of risk zone delineation and engineering practicality, and fundamentally reduces ineffective investment and risk omissions in water control projects.
[0095] The core objective is to transform the delineated risk areas into standardized results that can be directly connected to underground construction without secondary conversion. This solves the problem that traditional simulation results can only be viewed but not used. Traditional simulation results are mostly maps and text descriptions in reports. On-site construction personnel need to reconvert coordinates and recalculate construction parameters. Errors are very likely to occur during the secondary conversion process, and it cannot adapt to the dynamic progress rhythm of the mining face.
[0096] The target area coordinate set output in this stage is not a simple coordinate list, but a standardized deliverable package that fully conforms to coal mining industry standards and is directly compatible with on-site construction. Its core characteristics are as follows: Standardized coordinate system and format: The output coordinates adopt a spatial coordinate system that is completely unified with the mining area, including complete three-dimensional boundary information such as the type of each target area (water-rich anomaly zone / potential water-conducting channel), risk level, coordinates of plane boundary inflection points, and elevation of the roof and floor. It can be directly imported into the coal mine underground measurement system and drilling construction control system without any secondary conversion; Supporting engineering parameters: For different types of target areas, the core design parameters of the corresponding water control engineering are output simultaneously: For water-rich anomaly zones, the design coordinates, final hole depth, and final hole stratum of the advanced drilling verification hole are output; For potential water-conducting channels, the target area range, grouting hole spacing and depth are output, directly compatible with mainstream water control engineering such as underground advanced drilling, floor grouting reinforcement, and working face pre-grouting; Dynamic updates to adapt to mining rhythm: The target area coordinate set can be updated in real time with the updates of preceding fiber optic monitoring data, hydrological parameter optimization, and flow field simulation results. For example, with each cycle of advancement of the working face, the system can complete a full-process iteration and synchronously update the target area coordinate set, completely solving the industry pain points of traditional static results lagging behind the mining progress and being unable to cope with dynamic water hazard risks.
[0097] In this embodiment, identifying and delineating water-rich anomaly zones and potential water-conducting channels involves: calculating the spatial variability of the hydraulic gradient and permeability coefficient tensor in the high-precision aquifer level and flow field distribution prediction results; calculating the flow dominance factor for each grid cell based on the hydraulic gradient and permeability coefficient tensor, which characterizes the likelihood of water flow preferentially passing through; performing three-dimensional spatial clustering analysis on the flow dominance factor to identify continuous regions with high flow dominance factors as potential water-conducting channels; and identifying regions with hydraulic gradients below a preset threshold as water-rich anomaly zones.
[0098] It should be noted that: 3D full tensor hydraulic gradient calculation: Traditional methods only calculate the hydraulic gradient in the x and y directions of the plane, completely ignoring the gradient change in the vertical z direction. However, over 90% of underground water inrushes in coal mines originate from vertical water-guiding channels (such as floor faults connecting deep Ordovician limestone aquifers, or roof fracture zones connecting overlying aquifers). The vertical hydraulic gradient is a core indicator for determining the driving force of vertical water flow. This step, based on 3D flow field head data, calculates the partial derivatives of the head along the x, y, and z directions for each grid cell, generating a 3D hydraulic gradient vector to fully capture the driving force characteristics of planar runoff and vertical overflow. Anisotropic permeability coefficient tensor adaptation: The permeability coefficient used in this step is the 3D anisotropic permeability coefficient tensor obtained from previous optimization (containing permeability coefficients along the three principal axes Kx, Ky, and Kz), rather than the scalar equivalent permeability coefficient of traditional methods. The difference between the bedding and vertical permeability coefficients of layered rock strata in coalfields can be 2-4 orders of magnitude. Scalar parameters will completely smooth out the anisotropic seepage characteristics of the rock strata. Tensor calculations in this step can accurately restore the original laws of bedding runoff and vertical seepage resistance of layered aquifers, as well as the anisotropic abrupt changes after the destruction of tectonic zones.
[0099] By using geostatistical variability functions, the spatial abrupt changes in the hydraulic gradient and permeability coefficient tensor around each grid cell are quantified. This fundamentally addresses the weakness of traditional methods that cannot distinguish between normal lithological changes and anomalous structures when interpreting parameters at a single point: normal lithological facies transitions result in continuous and gradual parameter changes with low spatial variability; however, water-conducting structures such as faults, collapse columns, and fracture zones exhibit abrupt changes in parameters, resulting in extremely high spatial variability. Quantifying spatial variability allows for the early filtering out of false anomalies caused by normal lithological fluctuations, pinpointing the true parameter anomaly areas for subsequent risk identification and significantly reducing the misjudgment rate in subsequent calculations.
[0100] Traditional methods for identifying water channels in coalfields suffer from two major pitfalls: First, relying solely on high permeability coefficients as the criterion. However, high-permeability areas without hydraulic gradients and water flow are not water channels (e.g., closed high-permeability sandstone lenses), and traditional methods misclassify them as water-risk zones. Second, relying solely on high hydraulic gradients as the criterion. However, high-hydraulic-gradient areas with extremely low permeability coefficients are also not water channels (e.g., stress concentration zones in aquitards), and traditional methods produce numerous invalid interpretations.
[0101] The two approaches essentially separate the two core elements of groundwater seepage: permeability (the flow capacity of the medium itself) and driving capacity (the water flow dynamics determined by the head difference). The water flow advantage factor constructed in this step perfectly solves this core deficiency.
[0102] The flow dominance factor is an original composite quantitative index for coalfield seepage characteristics. Its core is tensor coupling calculation based on Darcy's law, and the formula is: flow dominance factor = |K·J|. Where K is the three-dimensional permeability coefficient tensor and J is the three-dimensional hydraulic gradient vector, the dot product of the two is the Darcy seepage velocity vector of the grid cell, and its magnitude is the flow dominance factor.
[0103] This indicator couples the flow capacity of the medium with the driving force of the water flow. The higher the value, the greater the groundwater flow rate through the unit per unit time, and the higher the probability that the water flow will preferentially pass through the unit, meaning it is more likely to become part of the water-conducting channel. Addressing the core risk of vertical water inrush in coal mines, amplified weights are specifically set for the vertical seepage component, strengthening the contribution ratio of the vertical water flow advantage factor. This can accurately capture hidden vertical water-conducting channels that traditional planar indicators cannot identify, perfectly meeting the core needs of preventing and controlling floor water inrush and roof water hazards in coal mines.
[0104] It transforms the previously qualitative and vague concept of water diversion channels in the field of coalfield water control into standardized numerical indicators that can be quantified, calculated, and compared for each grid unit, completely eliminating the subjectivity of traditional manual interpretation and providing a unique and unified quantitative benchmark for the accurate identification of subsequent water diversion channels.
[0105] Traditional methods for identifying water-conducting channels generally employ two-dimensional planar isoline delineation, which suffers from two major drawbacks: first, they cannot identify concealed water-conducting channels with vertical extension and complex spatial morphology (such as stepped fault combinations, concealed collapse columns with vertical water-conducting zones, and network-like fracture development zones), resulting in an extremely high false negative rate; second, they cannot distinguish between isolated high-value false anomalies and continuous water-conducting channels, leading to a high false positive rate. This approach addresses these two problems at their root through three-dimensional spatial clustering.
[0106] The DBSCAN density clustering algorithm, rather than the more common K-means clustering, is specifically chosen for its adaptation to the spatial characteristics of water-conducting channels. DBSCAN does not require pre-setting the number of clusters and can automatically identify continuous high-density regions of arbitrary shapes, perfectly adapting to the complex spatial morphology of water-conducting channels (fault zones are strip-shaped, collapse columns are columnar, and fracture zones are network-like). K-means, on the other hand, can only identify spherical clusters and is completely unsuitable for the irregular shapes of water-conducting channels. The clustering process employs two constraint rules: ① Numerical constraint: only grid cells with a flow dominance factor higher than a preset threshold are clustered (the threshold is taken as the 95th percentile or higher of the background value of the aquifer in the mining area, locking in cells with high flow capacity across the entire region); ② Spatial constraint: only continuously connected high-value cells in three-dimensional space are clustered into the same cluster. Isolated high-value cells (such as local high-permeability lenses) do not meet the spatial continuity requirement and will not be identified as water-conducting channels, completely filtering out false anomalies.
[0107] The clustering output is not a two-dimensional planar range, but a complete three-dimensional geological model of the water diversion channel. It accurately includes core engineering parameters such as the spatial orientation of the channel, vertical extension depth, connectivity range, spatial distance from the mining face / tunnel, and the aquifer it connects to. It can directly determine whether the channel poses a water inrush threat and the threat level, providing accurate spatial boundaries for subsequent water control engineering design.
[0108] A fatal flaw is commonly found in traditional methods of identifying water-rich anomalies in coalfields: directly equating high-permeability areas with water-rich anomalies. However, from a hydrogeological perspective, high-permeability areas, even if they are strong flow zones with large hydraulic gradients and continuously flowing groundwater, have limited storage capacity and are not the core source of mine flooding. In contrast, water-rich anomalies are essentially groundwater catchment / stagnant zones within aquifers with large storage space, high head pressure, and slow flow. Once exposed by mining operations, they instantly release large amounts of statically stored water, making them the core source of mine flooding. This section corrects this industry misconception by defining a low hydraulic gradient threshold.
[0109] The identification of water-rich anomaly zones in this stage adopts a dual-control rule of main indicators + auxiliary constraints to avoid misjudgment by a single indicator: The core main indicator is that the hydraulic gradient is lower than the preset threshold. The threshold is set based on the background hydraulic gradient of the aquifer in the mining area. Usually, 1 / 5 of the background hydraulic gradient of the area is taken as the critical value. Areas below this threshold are areas with slow water flow and strong water storage capacity. Auxiliary constraints are: the permeability coefficient and porosity are simultaneously superimposed to ensure that the identified area not only has slow water flow, but also has sufficient water storage space and flow capacity, and excludes invalid stagnant areas with low permeability and low porosity.
[0110] The occurrence of coal mine water inrush accidents requires the simultaneous fulfillment of two core conditions: sufficient water source (abnormally water-rich area) and a conductive channel (potential water-conducting channel). This solution uses two independent quantitative systems to accurately identify the water source and channel of the water inrush, comprehensively covering the two core elements of water hazard prevention and control. It completely solves the one-sided problems of traditional methods that only identify channels while ignoring water sources or only consider water abundance while ignoring water-conducting paths, thus achieving full-chain identification of water hazard risks.
[0111] In this embodiment, the parameter optimizer based on the physical information neural network further includes a step of dynamically adjusting the weights of the physical equation constraint terms in the loss function during the iterative optimization process. Specifically, this involves: calculating the average absolute error between the simulated water level and the measured water level in the current iteration step; adaptively adjusting the weight coefficients of the physical equation constraint terms according to the changing trend of the average absolute error; and increasing the weight coefficients of the physical equation constraint terms when the average absolute error tends to stabilize, so as to strengthen the constraint of physical laws.
[0112] It should be noted that the Mean Absolute Error (MAE) is chosen as the criterion for determining the convergence state of the fit, rather than the Root Mean Square Error (RMSE) used in the calibration phase of this scheme. This choice is specifically tailored to the characteristics of the PINN training process and coalfield monitoring data. The core reason is that RMSE is extremely sensitive to large deviations and outliers, and during training, it can oscillate violently due to local error fluctuations at individual monitoring points, failing to stably reflect the convergence trend of the network's overall fitting ability. In contrast, the water level data from observation wells in the coalfield inevitably contains occasional outliers caused by interference from underground equipment, temporary forced drainage, and monitoring malfunctions. MAE is not sensitive to outliers and can stably and accurately reflect the overall convergence state of the network's fit to the full set of measured data. MAE directly represents the average absolute deviation between simulated and measured values, with intuitive physical meaning. It directly corresponds to the absolute error requirements for water level prediction in coal mine water control, accurately pinpointing the bottleneck nodes of the network's fitting ability and providing a clear and reliable trigger basis for subsequent weight adjustments.
[0113] The error calculation in this step is strictly real-time and global, perfectly suited to the iterative training rhythm of PINN: Calculation frequency: In each iteration of network training, a full MAE calculation is performed synchronously, rather than fixed-interval sampling calculations, ensuring that weight adjustments can dynamically follow the training state in real time without lag. Calculation scope: Covering full-time measured water level data from all observation wells, simulated water level data from the benchmark model, and constraint data corresponding to the dynamic parameter field retrieved from fiber optics, rather than using only single-point data from a few observation wells, ensuring that the judgment results reflect the network's global fitting ability, rather than local fitting effects.
[0114] In traditional PINN applications, the weights of each component of the loss function are fixed values preset before training. In coalfield hydrological scenarios, this presents two intractable engineering deadlocks: **Early-Training Deadlock:** If the weights of the physical equation constraint terms are preset too high, the network will be locked by strong physical rules in the early stages of training, unable to quickly fit the measured data, resulting in extremely slow convergence or even complete non-convergence. If the weights are preset too low, the network will quickly overfit the noise and local biases of the measured data, getting stuck in local optima and outputting hydrological parameters without physical meaning, which cannot be corrected later. **Late-Training Deadlock:** When the data fitting error stabilizes, the fixed low weights cannot enhance the physical consistency of areas without monitoring data, leading to good fitting around observation wells and complete distortion of the overall large-scale prediction, completely lacking the generalization prediction ability necessary for coal mine water control. Setting high weights initially, however, leads to the early-stage non-convergence deadlock. Meanwhile, the core pain point in the coalfield scenario is that the number of monitoring wells is extremely small and their spatial distribution is extremely uneven (concentrated in mining areas, with no monitoring data in remote areas). Fixed weights are completely unable to meet the needs of having data in some areas and requiring constraints across the entire domain. This is also the core starting point for the design of this section.
[0115] The weight adjustments in this stage are not general linear increases or decreases, but are designed entirely around the core needs of coalfield hydrological simulation. They involve phased, boundary-based, and trigger-based custom adjustment rules, which are divided into three core stages: In the initial training phase (rapid MAE descent phase), the criterion is: if the relative decrease in MAE exceeds 5% within several consecutive iterations, it is considered a phase of rapid convergence of the fitting error. The adjustment strategy is to set the weights of the physical equation constraint terms to low initial values (typically 0.1-0.5), retaining only basic physical constraints to avoid strong physical rules limiting the network's fitting ability. The core objective of this phase is to enable the network to quickly inherit prior information from the preceding benchmark model and the dynamic parameter field retrieved from the fiber, anchoring it to a reasonable parameter space that conforms to the measured data, thus solving the problem of convergence difficulties in the early stages of traditional PINN training.
[0116] Mid-training (MAE decline slowing phase): The criterion is that if the relative decrease in MAE is between 0.1% and 5% over several consecutive iterations, it is considered a slowing-down phase of convergence. Adjustment strategy: A linear incremental strategy is adopted, gradually increasing the weight coefficients of the physical equation constraint terms by 5%-10% every 10 iterations, balancing data fitting accuracy with consistency with physical laws. The core objective at this stage is to avoid the network overfitting to noise in local monitoring data. While fitting the measured data, the global physical constraints are gradually strengthened to prevent the network from getting trapped in local optima.
[0117] In the later stages of training (when MAE stabilizes), the judgment rule is: if the relative decrease in MAE is less than 0.1% over 30 consecutive iterations, the fitting error is considered to have stabilized, and the data fitting has reached a bottleneck. The adjustment strategy is to significantly increase the weight coefficients of the physical equation constraint terms, typically to 5-10 times the initial value (not exceeding a preset weight limit), forcibly strengthening the global constraints of the physical laws governing groundwater flow. The core objective of this stage is to correct local deviations from hydrological laws that appeared in the earlier fitting, especially in blank areas without monitoring data. This forces the network's output of the head field and hydrological parameter field to satisfy the groundwater flow control equations across the entire domain, completely resolving the core pain points of traditional PINN's prediction distortion and poor generalization in areas without data.
[0118] To avoid extreme situations in weight adjustment, strict upper and lower limits are set for weights in this step: the lower limit of weights is no less than 0.1 to ensure that physical constraints are effective throughout the entire training process and that there will never be a black box fitting driven by pure data; the upper limit of weights is no more than 10 to avoid the data fitting accuracy from being reduced due to excessively high physical weights, and to always maintain the core bottom line that the accuracy of water level prediction meets the requirements of coal mine engineering.
[0119] This section is not a general algorithm optimization, but an original adaptation specifically for coalfield hydrological engineering scenarios. Its core innovative value lies in three aspects: First, it completely resolves the inherent contradiction between data fitting speed and physical consistency in PINN training, achieving a dynamic balance between the two. This allows PINN to stably and efficiently converge to the globally optimal hydrological parameter set, rather than a locally optimal solution, even in extreme scenarios with sparse monitoring data and highly heterogeneous aquifers in coalfields. Second, it achieves fully automated engineering implementation of PINN training, completely solving the industry pain point of traditional fixed-weight PINN requiring repeated manual trial and error parameter tuning over a period of several weeks. It eliminates the need for repeated debugging by hydrogeologists and algorithm engineers, significantly lowering the application threshold of PINN in coal mines. Third, it specifically addresses the pain point of uneven distribution of coalfield monitoring data. By strengthening the physical constraints across the entire domain in the later stages, the prediction results in remote areas without monitoring data and structurally complex areas strictly follow the groundwater flow patterns, ensuring the reliability of the overall flow field prediction and fully meeting the core needs of comprehensive risk prevention and control in coal mine water management.
[0120] Example 2, Figure 2 This invention presents a system for optimizing and numerically simulating hydrological parameters of coalfield aquifers, comprising: a data acquisition and processing module for acquiring and processing drilling monitoring data, historical hydrogeological data, and temperature and strain time-series data collected by a distributed fiber optic sensing system; a model construction and initialization module for constructing an initial aquifer structure model, an initial hydrological parameter field, and a baseline numerical simulation model; a parameter dynamic inversion module for inverting a dynamic hydrological parameter field based on temperature and strain time-series data; a physical information neural network optimization module, which incorporates a parameter optimizer based on a physical information neural network to generate an optimized set of hydrological parameters; a numerical simulation and prediction module for performing multi-scale simulation calculations based on the optimized set of hydrological parameters and outputting prediction results; and a visualization and target area delineation module for three-dimensional visualization and identification of water-rich anomaly areas and potential water-conducting channels.
[0121] The above formulas are all dimensionless calculations. The formulas are derived from software simulations based on a large amount of collected data to obtain the most recent real-world results. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.
[0122] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented using software, the above embodiments can be implemented, in whole or in part, in the form of a computer program product.
[0123] Those skilled in the art will recognize that the modules and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.
[0124] In addition, the functional modules in the various embodiments of this application can be integrated into one processing module, or each module can exist physically separately, or two or more modules can be integrated into one module.
[0125] The above are merely specific embodiments of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.
Claims
1. A method for optimizing and numerically simulating hydrological parameters of coalfield aquifers, characterized in that, Includes the following steps: Step S1: Obtain drilling monitoring data and historical hydrogeological data from multiple boreholes within the target coalfield area, construct an initial aquifer structure model, and generate an initial hydrological parameter field based on the initial aquifer structure model; Step S2: Based on the initial hydrological parameter field, construct an initial groundwater numerical simulation model, and use preset observation well water level data to calibrate and verify the initial groundwater numerical simulation model to obtain a benchmark numerical simulation model. Step S3: Deploy a distributed optical fiber sensing system in the target coalfield area to acquire real-time temperature and strain time-series data from multiple spatial locations, and invert the dynamic hydrological parameter field based on the temperature and strain time-series data; the inversion of the dynamic hydrological parameter field specifically involves: preprocessing the temperature and strain time-series data by denoising and normalization to extract its spatiotemporal variation characteristics; constructing an inversion model based on the heat conduction equation and the thermo-solid coupling constitutive relation, using the preprocessed temperature and strain time-series data as input to the inversion model, and using the aquifer permeability coefficient and porosity as the parameters to be inverted; using a Bayesian inversion framework combined with the Markov chain Monte Carlo method to solve the inversion model and obtain the posterior probability distribution of the parameters to be inverted; extracting the mean of the posterior probability distribution as the dynamic hydrological parameter field; Step S4: Construct a parameter optimizer based on a physical information neural network. Take the simulated water level of the benchmark numerical simulation model, the dynamic hydrological parameter field, and the water level data of the observation well as inputs. Iteratively optimize the hydrological parameters by solving the loss function embedded in the partial differential equation of groundwater flow, and generate an optimized set of hydrological parameters. Step S5: Input the optimized hydrological parameter set into the benchmark numerical simulation model, update the model parameters, perform multi-scale simulation calculations, and output high-precision prediction results of aquifer water level and flow field distribution.
2. The method for optimizing and numerically simulating hydrological parameters of coalfield aquifers according to claim 1, characterized in that, The construction of the initial aquifer structure model in step S1 specifically includes: The drilling monitoring data and historical hydrogeological data were cleaned and standardized to extract lithological stratification information and marker layer depths for each borehole. A sequential indicator simulation algorithm is adopted, using the lithological stratification information as hard data and the stratigraphic data interpreted by seismic exploration as soft data, to perform three-dimensional spatial interpolation and random simulation, generating multiple lithological distributions with equal probability. Calculate the average probability distribution of the multiple equally probable lithological distributions and use it as the initial aquifer structure model.
3. The method for optimizing and numerically simulating hydrological parameters of coalfield aquifers according to claim 2, characterized in that, The calibration and verification described in step S2 specifically include: The first time period data from the observed well water level data is extracted as calibration period data, and the second time period data is used as verification period data. The hydrological parameters of the initial groundwater numerical simulation model are automatically calibrated using a parallel particle swarm optimization algorithm with the goal of minimizing the root mean square error between the simulated and measured water levels in the calibration period data. The calibration model is verified using the data from the validation period. When the Nash efficiency coefficient during the validation period is greater than a preset threshold, the current model is determined as the benchmark numerical simulation model.
4. The method for optimizing and numerically simulating hydrological parameters of coalfield aquifers according to claim 3, characterized in that, The loss function obtained by solving the partial differential equation embedded in the groundwater flow in step S4 is specifically as follows: A deep neural network consisting of an input layer, multiple fully connected hidden layers, and an output layer is constructed. The input layer is used to receive spatial coordinates and time coordinates, and the output layer is used to output the predicted water head value. A total loss function is constructed, which includes data fitting terms, initial condition constraints, boundary condition constraints, and physical equation constraints. The physical equation constraints are constructed based on the groundwater flow continuity equation and Darcy's law, and are used to constrain the predicted hydraulic head value output by the deep neural network to satisfy the physical laws of groundwater flow. An adaptive moment estimation optimizer is used to iteratively update the network weights of the deep neural network with the goal of minimizing the total loss function. After training, the optimized hydrological parameters are extracted from the deep neural network.
5. The method for optimizing and numerically simulating hydrological parameters of coalfield aquifers according to claim 4, characterized in that, The multi-scale simulation calculation described in step S5 specifically includes: Construct a nested grid system, using locally finer grids in the well perimeter area and areas with complex geological structures, and coarser grids in the remaining areas; The optimized set of hydrological parameters is mapped to each grid cell of the nested grid system; Based on the nested grid system, the groundwater flow control equations are discretized using the finite volume method to construct a large-scale sparse linear equation set. The large-scale sparse linear equations are solved in parallel using the algebraic multigrid preprocessing conjugate gradient method, resulting in high-precision predictions of aquifer water level and flow field distribution.
6. The method for optimizing and numerically simulating hydrological parameters of coalfield aquifers according to claim 5, characterized in that, It also includes step S6: The high-precision aquifer water level and flow field distribution prediction results are rendered in three dimensions to generate a dynamic interactive display interface. Based on the dynamic interactive display interface, water-rich anomaly areas and potential water diversion channels are identified and delineated. The boundary information between the water-rich anomaly zone and the potential water-conducting channel is output as the target area coordinate set for coal mine water control work.
7. The method for optimizing and numerically simulating hydrological parameters of a coalfield aquifer according to claim 6, characterized in that, The identification and delineation of water-rich anomaly zones and potential water-conducting channels specifically involves: Calculate the spatial variability of the hydraulic gradient and permeability coefficient tensors in the high-precision aquifer level and flow field distribution prediction results; Based on the hydraulic gradient and permeability coefficient tensor, the flow dominance factor of each grid cell is calculated. The flow dominance factor is used to characterize the probability of water flow passing preferentially. Three-dimensional spatial clustering analysis was performed on the aforementioned flow dominance factors to identify continuous regions with high flow dominance factors as potential water-conducting channels; Areas with hydraulic gradients below a preset threshold are identified as water-rich anomaly zones.
8. The method for optimizing and numerically simulating hydrological parameters of coalfield aquifers according to claim 4, characterized in that, The parameter optimizer based on the physical information neural network further includes a step of dynamically adjusting the weights of the physical equation constraint terms in the loss function during the iterative optimization process, specifically: Calculate the mean absolute error between the simulated water level and the measured water level at the current iteration step; Based on the changing trend of the mean absolute error, the weight coefficients of the constraint terms of the physical equation are adaptively adjusted. When the mean absolute error tends to stabilize, the weight coefficients of the constraint terms of the physical equation are increased to strengthen the constraint of physical laws.
Citation Information
Patent Citations
Hydrological numerical simulation calculation method based on meshless calculation
CN116151152A
Pollutant migration key parameter acquisition method under site scale temperature-hydrodynamic force coupling influence
CN117408178A