Fractured stratum groundwater pollution source positioning method based on deep learning enhanced inversion

Abstract

The invention discloses a fractured stratum groundwater pollution source positioning method based on deep learning enhanced inversion, and relates to the field of pollution source positioning, and the method comprises the steps: constructing a space-time concentration response data set containing complex physical field characteristics based on hydrogeology and pollution source parameters, and training a deep learning substitution model according to the space-time concentration response data set; the complex groundwater solute transport process is efficiently approximated, so that the problem that traditional numerical simulation consumes too long time is solved. On the basis, an integrated inversion framework is constructed, time sequence concentration data observed on site is combined with a substitution model, and global optimization is carried out by utilizing a genetic algorithm. Not only can a nonlinear high-dimensional parameter space be processed, but also a space-time evolution rule in observation data can be fully excavated, finally, synchronous accurate inversion of pollution source three-dimensional coordinates, dynamic release history and key hydrogeological parameters is achieved, and the defects that a traditional method is low in recognition rate and poor in timeliness under complex geological conditions are effectively overcome.

Description

Technical Field

[0001] This application relates to the field of pollution source location, and more specifically, to a method and system for locating groundwater pollution sources in fractured strata based on deep learning-enhanced inversion. Background Technology

[0002] In fractured aquifers, groundwater flow and pollutant migration are strongly controlled by geological structures. The three-dimensional fracture network determines the migration path, velocity, and spatial distribution of pollutants, making pollution sources more concealed and extremely difficult to identify. Therefore, in order to effectively assess pollution risks and formulate remediation plans, it is particularly urgent to develop an efficient groundwater pollution source inversion and location scheme that can adapt to complex geological conditions.

[0003] Existing methods for identifying groundwater pollution sources mainly include geophysical exploration, geochemical analysis, and simulation optimization. Among these, the simulation optimization framework is favored because it can comprehensively invert the location and release history of pollution sources; however, its high computational cost limits its application in complex three-dimensional fractured aquifers. To address efficiency issues, existing techniques often introduce alternative models (such as support vector machines and conventional neural networks) to approximate the numerical simulation process. However, current inversion schemes based on alternative models still have significant drawbacks: First, existing studies are mostly based on the assumption of aquifer homogeneity, neglecting the significant heterogeneity of fractured aquifers, multi-scale coupling, and exchange interactions between fracture and matrix domains. This results in models failing to accurately reflect the preferred pathway migration behavior of pollutants in the fracture network, leading to low reliability of the inversion results. Second, despite the introduction of alternative models, existing deep learning applications often fail to fully extract time-series information from observational data. Many methods rely solely on static input-output mapping or time-averaged concentration data, failing to explicitly incorporate complete time-series changes into the inversion process. This masks key dynamic features reflecting the release history of pollution sources, making it difficult to accurately reconstruct time-varying release histories.

[0004] Therefore, existing technologies urgently need an inversion scheme that can integrate complex discrete fracture network modeling with high-precision temporal deep learning feature extraction. Summary of the Invention

[0005] To address the aforementioned problems in the existing technology, according to one aspect of this application, a method for locating groundwater pollution sources in fractured strata based on deep learning-enhanced inversion is provided, comprising: a training phase and a location phase; The training phase includes: constructing a spatiotemporal concentration response dataset based on a set of hydrogeological parameters and a set of pollution source parameters; and training and validating the alternative model based on the spatiotemporal concentration response dataset to obtain the trained alternative model. The localization phase includes: acquiring field observation concentration data; constructing an integrated inversion framework based on the field observation concentration data and a trained alternative model, the integrated inversion framework including decision variables, objective function, and constraints; performing global optimization and parameter inversion on the integrated inversion framework based on a genetic algorithm to obtain the optimal parameter solution; and performing comprehensive interpretation of pollution source characteristics on the optimal parameter solution to obtain a pollution source characteristic report, the pollution source characteristic report including the three-dimensional coordinates of the pollution source, time series data representing the release history, and estimated values ​​of various hydrogeological parameters.

[0006] According to another aspect of this application, a groundwater pollution source localization system based on deep learning-enhanced inversion is provided, which includes: a training module and a localization module; The training module includes: a spatiotemporal concentration response data construction unit, used to construct a spatiotemporal concentration response dataset based on a set of hydrogeological parameters and a set of pollution source parameters; and an alternative model training unit, used to train and validate an alternative model based on the spatiotemporal concentration response dataset to obtain a trained alternative model. The positioning module includes: an observation concentration data acquisition unit for acquiring on-site observation concentration data; an integrated inversion framework construction unit for constructing an integrated inversion framework based on on-site observation concentration data and a trained alternative model, the integrated inversion framework including decision variables, objective function, and constraints; an optimal parameter solution analysis unit for performing global optimization and parameter inversion on the integrated inversion framework based on a genetic algorithm to obtain the optimal parameter solution; and a pollution source feature report generation unit for performing comprehensive interpretation of pollution source features on the optimal parameter solution to obtain a pollution source feature report, the pollution source feature report including the three-dimensional coordinates of the pollution source, time series data representing the release history, and estimated values ​​of various hydrogeological parameters.

[0007] Compared with existing technologies, this application provides a method and system for locating groundwater pollution sources in fractured strata based on deep learning-enhanced inversion. Addressing the technical problems of high computational cost of numerical simulation, neglect of the heterogeneity of fractured aquifers, and low inversion accuracy due to insufficient utilization of time-series data in existing groundwater pollution source identification methods, this approach first constructs a spatiotemporal concentration response dataset containing complex physical field characteristics based on hydrogeological and pollution source parameters. A deep learning alternative model is then trained based on this dataset to efficiently approximate the complex groundwater solute transport process, thus solving the problem of excessively long processing times in traditional numerical simulations. Furthermore, an integrated inversion framework is constructed, combining time-series concentration data from field observations with the alternative model, and using a genetic algorithm for global optimization. This not only handles nonlinear high-dimensional parameter spaces but also fully explores the spatiotemporal evolution patterns in the observation data, ultimately achieving simultaneous and accurate inversion of the three-dimensional coordinates of pollution sources, dynamic release history, and key hydrogeological parameters. This effectively overcomes the shortcomings of traditional methods, such as low identification rate and poor timeliness under complex geological conditions. Attached Figure Description

[0008] The above and other objects, features, and advantages of this application will become more apparent from the more detailed description of the embodiments of this application in conjunction with the accompanying drawings. The drawings are provided to further illustrate the embodiments of this application and form part of the specification. They are used together with the embodiments of this application to explain this application and do not constitute a limitation thereof. In the drawings, the same reference numerals generally represent the same components or steps.

[0009] Figure 1 This is a flowchart of a method for locating groundwater pollution sources in fractured strata based on deep learning-enhanced inversion, according to an embodiment of this application.

[0010] Figure 2 This is an architecture diagram of an alternative model in the deep learning-based augmented inversion method for locating groundwater pollution sources in fractured formations, according to an embodiment of this application.

[0011] Figure 3 This is a release history curve diagram in the deep learning-based augmented inversion groundwater pollution source localization method for fractured strata according to an embodiment of this application.

[0012] Figure 4 This is a block diagram of a groundwater pollution source localization system based on deep learning-enhanced inversion in fractured strata, according to an embodiment of this application. Detailed Implementation

[0013] Embodiments of this disclosure will now be described in more detail with reference to the accompanying drawings. While some embodiments of this disclosure are shown in the drawings, it should be understood that this disclosure can be implemented in various forms and should not be construed as limited to the embodiments set forth herein. Rather, these embodiments are provided to provide a more thorough and complete understanding of this disclosure. It should be understood that the accompanying drawings and embodiments of this disclosure are for illustrative purposes only and are not intended to limit the scope of protection of this disclosure.

[0014] To address the problems in the existing technology, this application proposes a method for locating groundwater pollution sources in fractured strata based on deep learning-enhanced inversion. Figure 1 This is a flowchart illustrating a method for locating groundwater pollution sources in fractured formations based on deep learning-enhanced inversion, according to an embodiment of this application. Figure 1As shown in the embodiment of this application, the method for locating groundwater pollution sources in fractured strata based on deep learning-enhanced inversion includes a training phase and a location phase. The training phase includes: step S1, constructing a spatiotemporal concentration response dataset based on a set of hydrogeological parameters and a set of pollution source parameters; step S2, training and validating an alternative model based on the spatiotemporal concentration response dataset to obtain a trained alternative model. The location phase includes: step S3, acquiring field observation concentration data; step S4, constructing an integrated inversion framework based on the field observation concentration data and the trained alternative model, the integrated inversion framework including decision variables, objective function, and constraints; step S5, performing global optimization and parameter inversion on the integrated inversion framework based on a genetic algorithm to obtain the optimal parameter solution; step S6, performing comprehensive interpretation of pollution source features on the optimal parameter solution to obtain a pollution source feature report, the pollution source feature report including the three-dimensional coordinates of the pollution source, time series data representing the release history, and estimated values ​​of various hydrogeological parameters.

[0015] It is worth mentioning that the core of this application lies in constructing a high-precision alternative model to solve the problem of high computational costs in traditional numerical simulations. Traditional methods rely on massive amounts of observational data for iterative calculations. When facing geological environments such as fractured aquifers, which are highly heterogeneous and have complex flow fields, the computational load increases exponentially, making it difficult to meet the dual requirements of real-time performance and accuracy. By introducing a training phase, the powerful nonlinear mapping capabilities of deep learning networks are utilized to establish a rapid prediction channel from hydrogeological parameters and pollution source characteristics to concentration response. The purpose of training is to allow the alternative model to fully learn the physical laws of pollutant migration in fractured networks, especially the exchange interactions and spatiotemporal evolution characteristics between fractures and matrix domains. This allows the model to replace tedious numerical simulations with extremely high speed in subsequent inversion processes, achieving efficient and accurate global optimization.

[0016] In step S1, a spatiotemporal concentration response dataset is constructed based on the hydrogeological parameter set and the pollution source parameter set. It should be understood that in fractured aquifers, groundwater flow and pollutant migration are strongly controlled by the geometric characteristics of the fracture network (such as density, dip angle, and strike) and matrix permeability, exhibiting significant anisotropy and multi-scale coupling characteristics. If simple mathematical statistics are directly used or model training relies solely on sparse field observation data, it often fails to cover the full picture of pollutant migration under complex geological conditions, resulting in poor model generalization ability and distorted inversion results. Therefore, constructing a spatiotemporal concentration response dataset based on the hydrogeological parameter set and the pollution source parameter set aims to generate large-scale samples containing rich physical field information through high-fidelity numerical simulations (such as hybrid discrete fracture network / equivalent porous media model DFN / EPM). This captures the dynamic migration trajectory of pollutants in three-dimensional space under different parameter combinations, providing high-quality training samples for deep learning models that cover the fracture preferential flow effect and time-varying release history characteristics, ensuring sufficient physical consistency and predictive reliability of the alternative model in subsequent inversions.

[0017] In one possible implementation, step S1, based on the hydrogeological parameter set and the pollution source parameter set, constructs a spatiotemporal concentration response dataset, including: step S11, randomly generating a three-dimensional fracture network and discretizing the model of the hydrogeological parameter set to obtain a discretized DFN / EPM model; step S12, solving the physical field control equations and generating a global concentration field based on the pollution source parameter set to obtain a global spatiotemporal concentration field; step S13, extracting monitoring well concentration time series data from the global spatiotemporal concentration field based on the virtual monitoring well location coordinate set to obtain a multi-well concentration time series; step S14, performing parameter-response pair correlation on the hydrogeological parameter set, the pollution source parameter set, and the multi-well concentration time series to obtain a spatiotemporal concentration response dataset.

[0018] Specifically, the implementation process of step S1 is as follows: In one possible implementation, step S11 involves randomly generating a three-dimensional fracture network and discretizing the model based on the hydrogeological parameter set to obtain a discretized DFN / EPM model, including: step S111, randomly generating a discrete fracture network geometric entity set based on the hydrogeological parameter set to obtain a discrete fracture network geometric entity set; step S112, constructing a hybrid geometric model and performing computational mesh partitioning on the discrete fracture network geometric entity set to obtain a discretized DFN / EPM model.

[0019] Step S111 is fundamental to constructing a high-fidelity groundwater numerical model; it is a digital reconstruction process that transforms statistical laws into concrete geometric entities. The input to this step is a set of hydrogeological parameters, a rigorously defined set of numerical values ​​used to describe the statistical characteristics of the underground fracture network. Before implementation, it is necessary to clarify the source and composition of this parameter set. The hydrogeological parameter set is obtained based on field geological exploration data, such as using borehole television imaging technology to count the number and occurrence of fractures, measuring the trace length distribution of fractures through outcrop geological surveys, and estimating the connectivity and aperture of fractures through water pressure tests. This parameter set specifically includes, but is not limited to: fracture density (e.g., 500 fractures / km³, determining the number of fractures per unit volume), the distribution pattern of fracture center point locations (e.g., Poisson process distribution), and statistical distribution parameters of fracture size (e.g., the minimum radius in a truncated power-law distribution). Maximum radius and attenuation factor ), and the statistical distribution parameters of fracture orientation (such as the Fisher constant in the Fisher distribution). (Mean dip angle and mean strike). First, initialization is performed. A structured empty list or dataset, named Discrete Fracture Network Geometric Entity Set, is created in computer memory. This dataset stores the geometric attributes of each fracture object to be generated, including but not limited to 3D coordinates, radius, and normal vector. Then, the program proceeds to determine the number of fractures. Based on the preset fracture density (e.g., 500 fractures / km³) in the hydrogeological parameter set and the volume of the study area (e.g., a model size of 1000m × 1000m × 100m, i.e., a volume of 0.1km³), the total number of fractures to be generated is calculated. In this example, the total number of cracks =500 × 0.1 = 50. This value determines the number of iterations for subsequent loop generation. Next, the core loop generation phase begins. The program starts a loop from 1 to... The loop generates an independent crack geometry in each iteration i. Within a single iteration, the coordinates of the center point are generated first. This process strictly follows the Poisson process assumption, meaning that the spatial distribution of cracks is completely random and uniform, with no spatial clustering effect. To achieve this, a linear mapping formula is used to map standard uniformly distributed random numbers to the physical model space. The specific implementation formula is as follows: In this formula, This represents the absolute coordinates of the center point of the i-th fracture in three-dimensional space. and The effective spatial boundary for crack generation is defined. In the case of global generation, these values ​​correspond to the boundary coordinates of the study area (e.g., The range is [0, 1000]. The range is [0, 1000]. The range is [0, 100]. , , These are three independent random variables, uniformly sampled from the closed interval [0,1] using a computer's pseudo-random number generator. This mathematical expression ensures that the generated fracture center points appear with equal probability anywhere within the study area, thus conforming to the spatial homogeneity assumption in geostatistics. For example, if a certain sampling yields... =0.5, =0.1, =0.9, then the coordinates of the center point of the crack will be calculated as (500, 100, 90) meters. After determining the location of the crack, the crack size is generated. Cracks in nature often exhibit a pattern of more small cracks and fewer large cracks, a distribution usually described by a power-law distribution. To avoid generating physically unrealistic infinitely large or infinitely small cracks, a truncated power-law distribution model is adopted. The size generation is implemented using an inverse transform sampling method, and its control formula is as follows: In this formula, It is the calculated equivalent radius of the current crack. and These represent the lower and upper physical cutoff limits for the fracture radius, respectively. These two parameters are directly derived from the aforementioned set of hydrogeological parameters. For example, they are set based on geological surveys. =30 meters =150 meters. It is the attenuation factor of the power-law distribution, determined by truncating the power-law distribution model. >1, which controls the steepness of the crack size frequency distribution curve. The larger the value, the more small the cracks will be, such as 2.6. It is a newly generated uniformly distributed random number within the interval [0,1]. The mathematical essence of this formula is the inverse function of the cumulative distribution function (CDF), which takes a uniformly distributed probability value as input. It can directly analyze random variables that conform to a specific power-law probability distribution. This step ensures that the generated fracture network conforms to fractal geometry at the scale, accurately reflecting the multi-scale heterogeneity of the subsurface medium. Subsequently, fracture orientation is generated. The spatial orientation (i.e., attitude) of the fractures determines the dominant channel direction of groundwater flow and is a key factor influencing pollutant migration paths. According to the definition of hydrogeological parameter set, fracture orientation follows a Fisher distribution. The Fisher distribution can be understood as a normal distribution defined on a three-dimensional sphere, used to describe the dispersion of a set of vectors around a certain average direction. Its probability density function is given by the following formula: In the formula, and These represent the polar angle (deviation angle) and azimuth angle of the fracture surface normal vector relative to the average principal direction, respectively. This is the Fisher constant, a core statistical parameter in the set of hydrogeological parameters. It is determined through the Fisher distribution model and quantifies the concentration of fracture orientation, such as 30. When As the value approaches infinity, all generated crack normal vectors will be strictly parallel to the average direction; when When the value approaches 0, the fracture orientation will exhibit a completely random isotropic distribution. During implementation, by randomly sampling the aforementioned probability density function and combining it with the average fracture dip angle (e.g., 10°) and average fracture orientation (e.g., 10°) given in the parameter set, the specific normal vector of the current i-th fracture is calculated. This step mathematically simulates the control effect of tectonic stress fields on the fracture direction of rock mass during geological history. Once the position, size, and orientation calculations for a single iteration are complete, the program will perform a solid encapsulation operation. The calculated center point coordinates will be... ,radius The determined spatial normal vector is encapsulated into a separate data object (e.g., a class instance or struct containing attribute fields). This object is then appended to the discrete fracture network geometric entity set initialized at the beginning of the step. At this point, a complete fracture generation loop ends, and the program determines whether the preset total number has been reached. If the target is not reached, the random seed is reset or the random function is called again to proceed to the next iteration; if the target has been reached, the loop terminates. Finally, the discrete fracture network geometric entity set is output. This is a set containing... A complete dataset of fracture objects, where each element is precisely defined by its geometric parameters. This dataset not only strictly matches the input set of hydrogeological parameters in terms of statistical properties, such as density, size distribution, and attitude distribution, but also constructs a three-dimensional skeletal structure in terms of geometric morphology.

[0020] Step S112 is a crucial bridge connecting geostatistical features and physical field numerical calculations. First, the model domain size needs to be defined. The model domain size refers to the geometric boundary of the entire simulation area in three-dimensional space. This parameter is directly derived from the hydrogeological parameter set mentioned in previous steps or the study area determined by field exploration. In this embodiment, the model domain size is defined as 1000m × 1000m × 100m (length × width × height). This size defines the physical spatial boundary of the groundwater flow field and solute transport, i.e., the matrix extent of the equivalent porous medium (EPM). This size is obtained based on geological survey data from the actual engineering site, identifying a rectangular or irregular polyhedral region containing the potential impact range of the pollution plume. Next, the matrix geometry is constructed. In the geometric modeling kernel of professional numerical simulation software (such as COMSOL Multiphysics) or a standalone three-dimensional computer-aided design (CAD) environment, a continuous three-dimensional entity representing the rock matrix is ​​created based on the aforementioned model domain size. In this example, the program generates a cuboid geometric object with coordinates x∈[0,1000], y∈[0,1000], and z∈[0,100]. This cuboid represents a porous matrix domain with relatively low permeability and high porosity, serving as the primary carrier for subsequent groundwater storage and slow seepage. Next, the program proceeds to embedding fracture entities. The program reads the discrete fracture network geometric entity set. This set contains... There are 50 (for example) independent crack objects, each carrying the coordinates of its center point. ,radius And the orientation information determined by the normal vector. Traversing the set, for each fracture entity, the modeling engine reconstructs its shape in 3D space based on its geometric parameters. The fracture is idealized as a two-dimensional disk or an extremely thin plate with a certain thickness. For example, for the i-th fracture in the set, the program generates a disk surface inside the aforementioned matrix cuboid based on its center point (500, 100, 90), the calculated radius 50m, and the plane equation determined by the inclination and azimuth angles. This process is repeated until all fractures in the set are transformed into visual geometric patches. Based on this, the step of forming a hybrid geometry model is implemented. This is the key operation that integrates independent matrix entities with fracture patches. The program calls the embedding or composition functions in Boolean operations to handle the topological relationship between the fracture surface and the matrix. A continuous hybrid geometry model refers to a unified geometry that has undergone topological repair, in which the fracture surface is recognized as an internal boundary within the matrix, rather than an independent floating object. This model accurately describes the spatial relationship between the fracture network and the matrix, ensuring that the fractures and matrix are geometrically seamlessly connected. This guarantees that during subsequent physics field solutions, the fluid and solute can exchange mass between the fractures and the matrix, i.e., the coupling term. (Physical basis). If there are intersections between fractures, this step will also automatically calculate the intersection line to form a connected fracture network skeleton. After completing the geometric construction, the core meshing operation is performed. The above continuous hybrid geometric model is imported into the mesh generation module. The mesh generation module is an algorithm component in numerical simulation software used to discretize a continuous geometric domain into a finite number of computational units. Its specific architecture includes a geometric resolver, a size control function generator, and a meshing algorithm engine (such as Delaunay triangulation or forward propagation). The processing aims to transform continuous partial differential equations into a discrete system of algebraic equations. In order to achieve a balance between computational accuracy and efficiency, non-uniform meshing rules and parameters need to be set. Since fractures are the dominant channels for groundwater flow, with fast flow velocities and large gradients, while the flow velocity in the matrix is ​​extremely slow, a multi-scale meshing strategy is required. In the embodiments of this application, the following parameters are set: an extremely high resolution mesh is forced inside the fracture surface, at the intersection of fractures, and near the interface between fractures and the matrix, and the minimum element size is set to 1.5m. This ensures accurate capture of the complex flow details within the fracture network and the solute exchange flux between the fracture and the matrix. Conversely, in the matrix interior far from the fractures, where the physical field changes more gently, a coarser mesh is used, with a maximum element size of 35m. This adaptive meshing strategy based on geometric features effectively controls the total number of meshes, avoiding excessive computation. After initiating an automatic mesh generation algorithm (such as the Delaunay algorithm), the algorithm first generates a two-dimensional triangular mesh on the fracture surface, and then uses these triangles as constraints to advance into the matrix to generate a three-dimensional tetrahedral mesh. The algorithm automatically smooths the transition of mesh size, ensuring that the element size gradient from the fine region of 1.5m to the coarse region of 35m meets the numerical stability requirements (e.g., the growth rate does not exceed 1.2). Finally, a discretized DFN / EPM model is output. This model is no longer a smooth geometry, but a complex data structure composed of millions of nodes and elements such as tetrahedrons. It contains the spatial coordinates of all nodes. The model includes topological connections between elements (i.e., which node constitutes which element) and attribute labels indicating different physical domains (fracture domains or matrix domains). For example, the model might contain 100,000 nodes and 500,000 tetrahedral elements, where elements labeled as fractures have high permeability and elements labeled as matrix have low permeability. This digital computational model, with its high-quality mesh generation, possesses all the geometric and topological conditions directly applicable to solving the governing equations of groundwater flow and solute transport using the finite element method (FEM) or finite volume method (FVM).

[0021] In one possible implementation, step S12, which involves solving the physical field control equations and generating a global concentration field based on the pollution source parameter set to obtain a global spatiotemporal concentration field, includes: step S121, solving the coupled groundwater flow field of the discretized DFN / EPM model based on the hydrogeological parameter set to obtain the global water pressure distribution and the velocity field of the fracture domain and matrix domain; step S122, solving the temporal pollutant migration problem based on the velocity field of the fracture domain and matrix domain, the hydrogeological parameter set, and the pollution source parameter set to obtain a global concentration snapshot sequence; and step S123, aggregating the global concentration snapshot sequence to obtain a global spatiotemporal concentration field.

[0022] Step S121 forms the physical basis for constructing the spatiotemporal concentration response dataset. First, parameter assignment and boundary condition setting are performed. At this stage, it is necessary to retrieve physical parameters closely related to the flow field from the hydrogeological parameter set. These parameters originate from field hydrogeological survey experiments, such as permeability coefficients obtained through pressure water tests, storage coefficients obtained through porous pumping tests, or porosity data obtained through laboratory core analysis. Specifically, this parameter set includes: matrix permeability. For example, 5×10 - 11 m 2 Storage coefficient of fractures and matrix ( , ), fluid viscosity Take the viscosity and fissure aperture of water at room temperature Such as 4mm, and roughness correction factor For example, 0.025. During the assignment process, the program iterates through every mesh element in the discretized model. For three-dimensional tetrahedral elements labeled as the matrix, the matrix permeability is directly assigned. and matrix storage coefficient For two-dimensional triangular elements labeled as fractures, their equivalent hydraulic properties must first be calculated. Based on a modified form of the cubic law, the equivalent permeability of the fracture is calculated using the following formula. : In the formula, This represents the ability of a fractured medium to allow fluid to pass through, measured in square meters. The average physical aperture of the fracture, i.e., porosity; This is a dimensionless roughness correction factor used to correct the flow velocity hindrance effect of fracture wall roughness. Taking the values ​​in this embodiment as an example, if the fracture aperture... =0.004m, roughness correction factor =0.025, then the calculated equivalent permeability of the fracture is... =0.004 2 / 12×0.025≈3.33×10 -8m 2 This calculation shows that the permeability of the fracture is much higher than that of the matrix, which is 3.33 × 10⁻⁶. -8 m 2 Much greater than 5×10 -11 m 2 This quantitatively verified the physical characteristics of fractures as a dominant channel for groundwater flow. After the calculation was completed, the... Values ​​are assigned to all fractured mesh elements. Next, hydraulic boundary conditions are applied to the external geometric boundaries of the model. A constant head boundary condition is applied to the left boundary (X=0 surface) of the model domain (1000m×1000m×100m), setting the hydraulic pressure to correspond to a water level height of 110.0m; a constant head boundary condition is applied to the right boundary (X=1000 surface), setting the hydraulic pressure to correspond to a water level height of 10.0m; the remaining boundaries (front, back, top, and bottom) are set as flux-free boundaries. This 100m macroscopic head difference, i.e., 110m-10m, constitutes the fundamental driving force for the groundwater flow from left to right within the entire model domain, generating a macroscopic hydraulic gradient of approximately 0.1. Subsequently, the core process of discretizing the coupled control equations is introduced. The finite element method (FEM) is used to transform the continuous partial differential equations describing groundwater flow into a discrete system of algebraic equations. Since the model includes two regions with vastly different physical properties—fractures and matrix—a coupled set of equations is constructed to achieve mass exchange between them through coupling terms. For the fractured domain, the governing equation for groundwater flow is: In the formula, The fracture storage coefficient characterizes the ability of fractures to release water. The hydraulic pressure scalar field to be solved; For time, The dynamic viscosity of water, The gradient operator along the tangential direction of the fracture surface indicates that the fluid flows within the two-dimensional fracture surface. For the source and sink terms, this represents the rate of water exchange between the matrix and the fissures. For the matrix, the governing equation for groundwater flow is: In the formula, The matrix storage coefficient, This is a standard three-dimensional gradient operator, indicating fluid flow within a three-dimensional matrix pore; the right side of the equation... With the fracture domain equation These are opposites and strictly follow the law of conservation of mass, meaning the amount of water flowing out of the matrix equals the amount of water flowing into the fissures. Coupling terms The value is determined by the pressure difference at the interface between the fracture and the matrix, expressed as: ,in The conduction coefficient ensures the hydraulic continuity of the two systems. During discretization, the Galerkin weighted residual method is used to integrate the above equations over a global grid. The infinite-dimensional pressure field is approximated as a linear combination of nodal pressures using shape functions, ultimately assembling into the form of… A large set of sparse matrix equations, in which This is the total stiffness matrix (which includes permeability information). This is the pressure vector for all nodes in the global domain. The load vectors are generated for the boundary conditions and source / sink terms. Finally, the equations are solved and the flow velocity is calculated. An efficient linear algebra solver (such as PARDISO or the MUMPS direct solver) is invoked to solve the above matrix equations. The output of the solution is the hydraulic pressure value at every grid node (e.g., hundreds of thousands of nodes) in the model. Mapping these discrete pressure values ​​back to geometric space yields the global water pressure distribution. This is a continuous scalar field that visually demonstrates the decreasing pressure from upstream (the pressure corresponding to a 110m head on the left) to downstream (the pressure corresponding to a 10m head on the right), while also reflecting the local pressure gradient distortion caused by the presence of fractures. Based on the obtained pressure field, Darcy's law is used for post-processing to calculate the velocity field. For each fracture element, its tangential pressure gradient is calculated. Multiplying this by the fracture conduction coefficient yields the flow velocity in the fracture domain. For each matrix element, calculate its three-dimensional pressure gradient. The matrix domain flow velocity was obtained. The final output consists of two independent vector fields: the velocity fields of the fracture domain and the matrix domain. Fracture domain velocity field It exhibits extremely high flow velocities, such as m / day, and its direction is strictly distributed along the fracture direction, forming a highway for pollutant migration; while the matrix domain flow velocity field The values ​​are relatively small, such as in the mm / day range, and the direction is controlled by the regional hydraulic gradient.

[0023] Step S122 is the core physical process simulation step in constructing the dataset. This step follows directly from the steady-state flow field data obtained in the previous steps, introducing a time dimension to simulate the dynamic diffusion and reaction processes of pollutants driven by the flow field. First, the input parameter set needs to be defined. The hydrogeological parameter set retrieved at this point includes not only the aforementioned flow field parameters but also parameters closely related to solute migration. These parameters are obtained through laboratory core testing or tracer experiments, specifically including: fracture porosity. For example, 0.7 and matrix porosity For example, 0.3; these two parameters determine the size of the space in which the medium contains the fluid; fluid diffusion coefficient. ,like and ,like The molecular diffusion capacity of solutes in fissure free water and matrix pore water are described, respectively; and the linear equilibrium partition coefficients describing adsorption are described. , and effective volume density of the medium , In addition, the mass transfer factor needs to be defined. This is used to quantify the mass transfer rate driven by concentration difference between the fracture and the matrix. Simultaneously, a set of pollution source parameters is acquired, defining the spatiotemporal characteristics of the pollution sources, and is set through random sampling during training data generation. For example, the three-dimensional spatial coordinates of the pollution sources are defined. Let x = 500m, y = 100m, z = 90m be a point in the left region of the model, and define a release flux function that varies with time. For example, a sinusoidal function with a period of approximately 1000 days and a peak value of 1 mol / s is used to simulate a non-steady leakage process. The implementation process begins with initialization and source term loading. The program first sets the initial concentration field of the entire model domain (including all fractures and matrix nodes) to zero, representing a background state before contamination occurs. Subsequently, based on the coordinates in the contamination source parameter set... The specific grid cell containing the coordinates is located in the discretized DFN / EPM model using a grid search algorithm. The flux is then released in a time series. This is loaded as a source-sink term into the governing equations of the unit. This means that at each computational time step, the unit releases a certain mass of pollutants into the surrounding environment. Next, retention and exchange parameters are calculated. When pollutants migrate in groundwater, they are often adsorbed by rock particles, resulting in a hysteresis effect. To quantify this effect, the retention factor is calculated based on a set of hydrogeological parameters. For fractured domains, the formula is used... Calculate the fracture retention factor For the matrix domain, using the formula Calculate the matrix retention factor These two dimensionless parameters, being greater than 1, directly reflect the proportion of the delay in pollutant migration velocity relative to groundwater flow velocity. For example, if... =2, meaning the average velocity of pollutants within the matrix is ​​only half that of the water flow. Following this, the core time-domain iterative solution phase begins. The time step is set. For example, with a simulation period of 1 day and a total simulation duration of 1000 days, the simulation is progressively advanced starting from t=0. Within each time step, the convection-dispersion equation solver (integrated into finite element software such as COMSOL) is invoked. For the fracture domain, the following governing equations are solved: The first term on the left-hand side of the equation describes the rate of change of concentration over time (taking into account the retention effect); the second term describes the rate of change of concentration due to flow rate. Driven convection and by diffusion coefficient Driven diffusion; on the right side The term represents the mass loss due to diffusion into the matrix domain. Simultaneously, for the matrix domain, the corresponding governing equations are solved: The right side here A positive value represents the mass replenishment obtained from the fracture domain. These two equations are solved by exchanging terms. To achieve strong coupling, where This is the mass transfer coefficient. This term indicates that when the fracture concentration... Greater than the matrix concentration At the same time, pollutants diffuse from the fissures into the matrix (the matrix acts as a storage source); conversely, when pollutants in the fissures are washed away, causing a decrease in concentration, pollutants stored in the matrix are released back into the fissures (the matrix acts as a secondary pollution source). The solver utilizes the velocity field obtained in the preceding step S121. and Given these conditions, the new concentration values ​​for all nodes in the global domain at the current time step are calculated using numerical integration. After each time step is solved, a concentration snapshot is stored. Concentration values ​​of all hundreds of thousands of nodes in the model (including fracture concentration) and matrix concentration The data is extracted and saved as a three-dimensional data matrix, called a concentration snapshot. As the time step progresses, for example from day 1 to day 1000, a series of such snapshots will be generated. The final output global concentration snapshot sequence is a data set arranged in chronological order, which fully records the entire process of the contamination plume from its source generation, rapid diffusion through the fracture network, infiltration and retention in the matrix, and decay over time.

[0024] Step S123 integrates discrete time slices into a unified high-dimensional data object, facilitating subsequent data mining and deep learning training. In implementation, a four-dimensional data structure is first constructed. A large four-dimensional array (or tensor) is allocated in computer memory, with its four dimensions defined as spatial coordinates. Spatial coordinates Spatial coordinates and time The size of the array is determined by the number of grid nodes and the number of time steps. Next, data filling is performed. The program iterates through the global concentration snapshot sequence. For each time point in the sequence... A snapshot, which contains spatial concentration distribution data, is mapped and populated into the corresponding four-dimensional array. Within the slice. In this way, originally scattered files or data blocks are integrated into a logically coherent whole. The resulting global spatiotemporal concentration field is a structured four-dimensional dataset, which allows for... The quadruple index allows for direct and fast lookup of pollutant concentration values ​​at any location within the model domain at any time. This data structure not only contains physical field information for the entire domain but also preserves complete temporal evolution characteristics.

[0025] Step S13 is a crucial dimensionality reduction process that transforms the global physical field data into sparse observation data. This requires obtaining a virtual monitoring well location coordinate set, a precisely defined set of three-dimensional coordinate points that simulates the layout of monitoring wells in actual engineering projects. The layout of the monitoring wells significantly impacts the reliability of the inversion results; a layout parallel to the water flow direction yields the highest inversion accuracy. Therefore, in this embodiment, the virtual monitoring well location coordinate set is set to include three key coordinate points, such as monitoring wells arranged parallel to the flow direction. , , If the model domain is 1000m×1000m×100m, and the pollution source is located at (500, 100, 90), then the coordinates of the monitoring well may be set as follows: Located upstream or near the source (400, 100, 90), Located in the middle range (600, 100, 90), Located downstream (800, 100, 90), this system captures the complete evolution of the pollution plume along the main stream direction. During implementation, spatial positioning is performed first. The program iterates through the aforementioned set of virtual monitoring well coordinates. For each specified monitoring well coordinate... The algorithm performs nearest neighbor search within the spatial grid of the global spatiotemporal concentration field. Because the discretized DFN / EPM model uses an unstructured grid (such as tetrahedral elements), the actual grid nodes may not fall precisely on the preset coordinates. For example, for a monitoring well with coordinates (600, 100, 90)... The algorithm calculates the Euclidean distance between the well and its surrounding grid nodes, locking the closest spatially located grid node (e.g., node ID #50234, coordinates 600.2, 99.8, 90.1) as the surrogate node for that monitoring well. Time series extraction is then performed. Once the spatial node representing the monitoring well is locked, the program slices and extracts data from that node along the time dimension t. From the initial simulation time t=0 to the final time t=1000 days, pollutant concentration values ​​for that node are read at fixed time steps (e.g., 1 day or 10 days). This process is like inserting a virtual probe into a four-dimensional data volume to extract a concentration curve that changes over time. For example, for a monitoring well... The extracted sequence may show: a concentration of 0 from t=0 to 200 days (the pollutant has not yet arrived); a non-zero reading starting at t=200 days (breakthrough moment); and a peak concentration of 20 mol / m³ at t=500 days. 3It then slowly descends. This curve is the penetration curve of the monitoring well. Next, sequence assembly and data structuring are performed. These will be respectively... , , The three extracted independent concentration time series are processed. If the total simulation duration is 1000 days and the output step size is 10 days, then each series contains 100 data points. The program integrates these three series into a unified dictionary or matrix structure to form a multi-well concentration time series. This object is a structured dataset whose keys are the monitoring well numbers and whose values ​​are the corresponding time series arrays.

[0026] Step S14 is the final synthesis step in constructing the deep learning training sample library. It rigorously pairs the causes (input parameters) that lead to the physical process with the effects (observed responses), thus providing data support for the supervised learning algorithm. First, sample construction is performed. The program integrates all the key input parameters used to drive this numerical simulation into a high-dimensional input feature vector. This feature vector contains 14 key variables: pollution source feature parameters (source coordinates) and leakage rate Characteristic values), and hydrogeological parameters (including matrix permeability). Matrix porosity Crack density fracture dip angle , fracture orientation , fracture aperture (etc.). For example, the input feature vector for a simulation might be [500, 100, 90, 1.0, 5 × 10⁻⁶]. -11 [0.3,500,10,10,0.004,...]. Simultaneously, the multi-well concentration time series is used as the corresponding output label. If a sliding window strategy is used to process time-series data, the output labels may be formatted as a concentration matrix within the corresponding time window. Then, association matching is performed. A data record object is created, and... and Bind them to form a complete parameter response pair. This pair of data represents the concentration response of a specific observation well under specific geological conditions and pollution source scenarios. This deterministic physical relationship is the foundation upon which deep learning alternative models can learn intrinsic patterns through training. To construct a statistically representative large-scale dataset, iterative generation and assembly operations are required. By writing automated scripts, the program controls the entire process (steps S11 to S14) to run repeatedly hundreds or thousands of times. At the beginning of each new loop, the program randomly samples according to a preset parameter distribution range to generate a new set of input parameters. For example, fracture density is 400-600 / km².3 Uniform sampling was conducted between 8° and 12° of fracture dip angle, and the matrix permeability was between 4×10⁻⁶. -11 Up to 6×10 -11 m 2 The data is sampled according to a log-normal distribution, while the pollution source leakage rate is sampled according to a sinusoidal function distribution. Each newly sampled set of parameters triggers a completely new process of generating a 3D fracture network, solving the flow field, simulating solute transport, and extracting monitoring well data, thus generating a new independent sample. Finally, the dataset is completed. After the preset number of iterations (e.g., 1000 or 2000), the program will process all generated samples. The data is collected and stored to form the final spatiotemporal concentration response dataset. This dataset is a large-scale structured file (such as HDF5 or CSV format), with each row corresponding to an independent numerical simulation experiment. The dataset includes various geological and pollution scenarios ranging from homogeneous to highly heterogeneous and from simple to complex, and covers the control effect of fracture network geometric parameters (such as strike and dip angle) on pollutant migration paths.

[0027] In step S2, the alternative model is trained and validated based on the spatiotemporal concentration response dataset to obtain the trained alternative model. Specifically, the alternative model here includes a local temporal feature extraction module, a gated sequence encoding module, a long dependency modeling module, a key information weighting module, and a decoding module. Figure 2 This is an architecture diagram of an alternative model in the deep learning-based enhanced inversion method for locating groundwater pollution sources in fractured aquifers according to an embodiment of this application. Correspondingly, when dealing with the highly complex inversion problem of identifying groundwater pollution sources in three-dimensional fractured aquifers, although the hybrid discrete fractured network / equivalent porous media (DFN / EPM) model constructed in the preceding steps can provide high-fidelity physical field simulation, its computational process involves solving a large system of partial differential equations, resulting in a lengthy single run. Within the inversion framework, it typically requires tens of thousands of iterative searches combined with global optimization algorithms. If the numerical model is directly called in each iteration, the computational cost will increase exponentially, making it impossible to meet the timeliness requirements of emergency response in practical engineering applications. Furthermore, the high heterogeneity of fractured media leads to complex nonlinear spatiotemporal evolution characteristics in pollutant migration, making it difficult for traditional simple regression models to capture this multi-scale dynamic dependency. Step S2 is used to build a high-precision alternative model based on a deep learning architecture. By using the spatiotemporal concentration response dataset generated in the previous steps for supervised training, the model can learn and internalize the complex nonlinear mapping law between hydrogeological parameters, pollution source characteristics and observation well concentration response, thereby generating a prediction tool that can replace expensive numerical simulation calculations at millisecond speeds.

[0028] In one possible implementation, step S2 is carried out as follows: The implementation process mainly relies on a specially designed deep learning hybrid architecture - CNN-GRU-LSTM-Attention model, which aims to capture local mutation features, long-term dependencies and weight information at key moments at the same time.

[0029] First, data preprocessing and partitioning are performed. The program reads the spatiotemporal concentration response dataset stored on the hard drive and randomly partitions it into training and test sets in a 7:3 ratio. For example, if the total number of samples in the dataset is 2000, 1400 samples are randomly selected for updating model parameters (training), and the remaining 600 samples are used to evaluate the model's generalization ability (testing). To accelerate the convergence of the gradient descent algorithm and eliminate the influence between parameters of different dimensions, the input features (parameters) and output labels (concentration sequences) are normalized. Input features include the location of pollution sources. Leakage rate and hydrogeological parameters (such as matrix permeability) The system includes 14 variables, such as crack density. The max-min normalization method is used to linearly map these data to the [0,1] interval. For time series data, a sliding window technique is employed to transform static parameter features into a dynamic sample format suitable for time series prediction. The sliding window length is set to 5, meaning the model uses information from the current and four past time steps to make predictions at each time step, thus generating a sample with the shape of... The three-dimensional input tensor.

[0030] Subsequently, the data enters the local temporal feature extraction module. This module consists of a one-dimensional convolutional neural network (1D-CNN). The processed input tensor is fed into the CNN layer, which uses a set of learnable convolutional kernels to slide across the time axis. The convolution operation aims to extract local patterns in the input sequence, particularly high-frequency features that reflect abrupt changes, inflection points, or oscillations in pollutant concentration breakthrough curves. For example, when the pollution plume front just reaches the monitoring well, the concentration rises sharply from zero, and the CNN can effectively identify this local change in gradient. After convolution operations and processing with nonlinear activation functions (such as ReLU), a set of local temporal feature vectors containing rich local details is output.

[0031] Next, the feature data flows to the gated sequence encoding module, namely the gated recurrent unit (GRU). GRU is a variant of recurrent neural networks (RNNs) specifically designed to solve the vanishing gradient problem in long sequence training. In this embodiment, the GRU receives local features from the CNN output and uses its internal gating mechanism to filter historical information. The core of the GRU lies in two gates: the update gate... and reset door Update Gate The calculation formula is: Reset the door The calculation formula is: in, It is the input at the current moment (features extracted by CNN). It is the hidden state from the previous moment. The Sigmoid activation function compresses the output to [0,1]. and This is the weight matrix that the model needs to learn during training. (Reset gate) How much of the past information was determined It needs to be forgotten. If... If the value approaches 0, the previous state is ignored, which helps the model quickly adapt to new patterns when faced with non-stationary hydrogeological changes. The candidate hidden states are calculated under the reset gate. : Ultimately, the hidden state By the update gate Adjusting historical state With candidate state The mixing ratio yields: This mechanism enables GRU to efficiently retain historical concentration trend information that is useful for current predictions while simplifying the network structure, and output gated features.

[0032] Subsequently, the data is transmitted to the long-term dependency modeling module, namely the Long Short-Term Memory (LSTM) network. Although the GRU already possesses a certain memory capacity, for processes like groundwater pollution that last for years or even decades, such as the 1000-day simulation period in this case, a more robust long-term memory mechanism is needed. LSTM achieves this by introducing independent memory unit states. LSTM contains three gates: the input gate... Forgotten Gate and output gate Forgotten Gate Determine from the cell state Which outdated information (such as early, minor concentration fluctuations) is discarded during the input gate? The cell state determines which new information (such as the latest changes in pollutant source emission fluxes) will be updated. The update combines forgetting and input operations, and the output gate control outputs how much information as the hidden state based on the current cell state. The corresponding weight matrix and bias vector are randomly initialized at the start of training using methods such as Xavier and are continuously adjusted during backpropagation. The LSTM module accurately captures the long-term hysteresis effects of pollutant migration in the fractured network (such as the adsorption and re-release of pollutants by the matrix), outputting long-term dependent sequences.

[0033] To further improve prediction accuracy, a key information weighting module, namely the Bahdanau attention mechanism, is introduced. Since information from different time steps in a long sequence contributes differently to the current prediction result (for example, the peak time of pollutant arrival at the monitoring well is far more important than the background value), the attention mechanism strengthens key information by dynamically assigning weights. First, the attention score is calculated. This measures the correlation between the current hidden state and the historical state. ,in, For splicing operations, for Activation function The scores are then normalized using the Softmax function to obtain the attention weights. Weight The value of is between 0 and 1, and the sum is 1. If a certain time t corresponds to the peak arrival time of the pollution plume, then the value of at that time is... This will significantly increase the complexity. Finally, the hidden states at all time steps are summed using weighted averages to generate an attention-weighted context vector. : This process enables the model to automatically focus on the critical moments that play a decisive role in concentration evolution, effectively solving the problem of information loss in long sequences.

[0034] Finally, decoding and prediction are performed. Context vector The input is fed into a fully connected layer for linear transformation, ultimately outputting a predicted value. This output is the predicted concentration time series, representing the model's estimate of future concentration changes in the three monitoring wells under given input parameters. The entire model training and validation process is driven by a model optimization algorithm. The loss function is set as the root mean square error (RMSE), used to quantify the difference between the predicted concentration and the actual concentration in the spatiotemporal concentration response dataset. ,in These are actual observations. These are the model's predicted values. The sample size is given. The Adam optimizer is used for parameter updates. The initial learning rate is set to 0.001, the batch size to 32, and the number of training epochs to 100. In each iteration, the gradient of the loss function with respect to all weights and biases in the model is calculated, and these parameters are updated using backpropagation to minimize the RMSE. To prevent overfitting, the learning rate is decayed to 30% of its original value every 25 epochs. After training, the trained alternative model is evaluated using a reserved test set. Metrics such as the coefficient of determination, root mean square error (RMSE), and mean absolute error (MAE) are calculated. The trained model with the final optimized weights is output as the trained alternative model. In step S3, on-site observation concentration data is acquired. It is understood that the objective of this application is to address the problem of highly concealed and difficult-to-identify pollution sources in three-dimensional fractured aquifers. Although the preceding training phase constructed a high-precision alternative model capable of rapidly predicting concentration responses through numerical simulation and deep learning, this merely establishes a tool for forward prediction. To truly solve the problem of unknown pollution sources in actual sites, it is necessary to incorporate real-world observation information into the framework and use known results to infer unknown causes. Therefore, introducing the positioning phase and executing step S3 to acquire on-site observation concentration data is to introduce real physical field information into the inversion system, providing a calibration benchmark and objective function calculation basis for subsequent global optimization, thereby driving the inversion algorithm to search for the pollution source characteristics and hydrogeological parameter solutions that best conform to objective facts in the parameter space.

[0035] In one possible implementation, step S3 is as follows: First, a network of monitoring wells is deployed in the actual contaminated fractured aquifer site based on hydrogeological survey results. To obtain the highest inversion reliability, the monitoring wells should be arranged parallel to the groundwater flow direction as much as possible. For example, in a 1000m × 1000m area, three monitoring wells (numbered P1, P4, P7) are arranged along the main flow direction, with their coordinates set as (400, 100, 90), (600, 100, 90), and (800, 100, 90), respectively. Using high-precision water quality sensors or periodic manual sampling and analysis, pollutant concentration data in each monitoring well are collected at a fixed frequency (e.g., once every 10 days) within a preset time period (e.g., 1000 days). The collected raw data needs to be cleaned and formatted to construct on-site observation concentration data. This is a structured two-dimensional array matrix. Its dimensions are M×T, where M=3 represents the number of monitoring wells and T=100 represents the number of time steps. For example, the elements in the matrix... This indicates the pollutant concentration value observed in monitoring well No. 4 on day 500, such as 20.5 mol / m³. 3 .

[0036] In step S4, an integrated inversion framework is constructed based on field observation concentration data and a trained alternative model. This framework includes decision variables, an objective function, and constraints. It should be understood that during the training phase, high-fidelity numerical simulations and advanced deep learning models establish a rapid prediction path from hydrogeological parameters to concentration response; however, this only constructs a forward predictor. When facing unknown pollution sources and complex hydrogeological conditions in actual sites, a reverse inference problem needs to be solved—that is, using limited observational data to infer the characteristics and geological parameters of potential pollution sources. Due to the high nonlinearity and ambiguity of groundwater systems, direct inversion is almost impossible. Therefore, introducing the integrated inversion framework and executing step S4 to construct decision variables, an objective function, and constraints transforms the inversion problem into a standard mathematical optimization problem, thereby achieving precise location and characterization of unknown pollution events.

[0037] In one possible implementation, step S4 is as follows: To implement this step, it is first necessary to define the mathematical form of the inversion optimization problem, including the vectorization of decision variables, the construction of the objective function, and the setting of constraints.

[0038] First, decision variables are vectorized. This process is based on an identified list of sensitive parameters. This list was determined during the preliminary sensitivity analysis phase and explicitly identifies unknown parameters that significantly influence pollutant migration processes and are to be retrieved. In this embodiment, the list of sensitive parameters includes 14 key variables: the three-dimensional spatial coordinates of the pollution source. Leakage rate parameters characterizing release history and key hydrogeological parameters (such as matrix permeability) Crack density fracture dip angle , fracture orientation Matrix porosity The program will process these 14 scalar parameters in a predetermined order (e.g., ...). , , , , (,...) are integrated and arranged into a single, 14-dimensional vector of decision variables. This vector In the algorithm, it represents a potential solution (individual), that is, a possible pollution scenario hypothesis.

[0039] Secondly, the objective function is formalized. Objective function The root mean square error (RMSE) is the sole criterion for evaluating the quality of the inversion results, and its design aims to quantify the difference between the model's predictions and the actual observations. This method uses RMSE as the objective function, and its mathematical expression is as follows: In the formula, It is the input decision variable vector; This is the total number of monitoring wells, which is 3 in this example (i.e., P1, P4, P7). It is the total number of observation points, for example, 100 time steps; The field observation concentration data obtained in step S3 represents the actual concentration value of the i-th monitoring well at time j. When the input parameter is At that time, the corresponding concentration value is predicted by calling a trained alternative model. The calculation logic of this function is: for each possible combination of parameters... The alternative model quickly generates a set of predicted concentration fields, and the inversion framework compares them point by point with the field measured data, calculates the mean of the sum of squared errors, and then takes the square root. The smaller the value, the better the current parameter combination. The closer it gets to the actual state of groundwater pollution.

[0040] Next, constraints are set. To prevent the optimization algorithm from searching in a physically meaningless parameter space (e.g., negative permeability or coordinates outside the study area), the decision variable vector needs to be based on prior knowledge boundaries. Each component in Upper and lower bounds are defined. Prior knowledge boundaries are derived from geological exploration reports or literature. The constraints set in this embodiment include: pollution source coordinates. rice, rice, meters; leakage rate mol / s; matrix permeability m²; fracture density Stripes / km³; Fract dip angle These constraints are in the form of inequalities. Fixed within the inversion framework, the genetic algorithm is forced to generate and mutate populations only within these physically plausible hypercubes.

[0041] After completing the above definitions, we proceed to the computational framework construction phase of the integrated proxy model and optimization engine. First, we encapsulate the fitness function `fitness_function`. This function serves as the interface connecting the deep learning model and the genetic algorithm. Internally, the data flow strictly follows these logics: A. Input data preprocessing: The function receives candidate solution vectors generated by the genetic algorithm. Because normalized data was used in the preceding steps of training the alternative model, therefore First, perform the same preprocessing. The program then calls the normalization parameters saved during the training phase, i.e., the maximum value of each feature. and minimum value ), using formula The parameters with physical dimensions are converted into dimensionless values ​​in the interval [0,1]. These normalized parameters are obtained through statistical analysis and persistently stored when the training set is constructed. B. Call the surrogate model for positive prediction: The normalized vectors... The input is fed into a pre-loaded, trained alternative model. This model utilizes its learned weights and biases to output predicted concentration time series within milliseconds. C. Calculate the objective function value: Using the RMSE formula above, calculate the predicted sequence. With global constants Error value between D. Fitness value transformation: Since the genetic algorithm seeks the maximum value of the fitness function, and the goal is to minimize the RMSE, a reciprocal transformation is used: ,in It is a very small positive number, such as 1e-6, used to prevent division by zero errors. This transformation ensures that solutions with smaller errors have higher fitness, thus having a greater chance of being preserved during evolution.

[0042] Finally, configure the genetic algorithm engine. Instantiate a genetic algorithm (GA) solver in the Python environment (e.g., using the DEAP or Geatpy library). Register the encapsulated `fitness_function` as the solver's evaluation core. Configure the constraint set as the boundary limits of the variables. Set the population size to 30, the crossover probability to 0.9 (using a single-point crossover strategy), the mutation probability to 0.2 (using polynomial mutation), and the maximum number of iterations to 500 generations. At this point, a complete ensemble inversion framework is built.

[0043] In step S5, the integrated inversion framework undergoes global optimization and parameter inversion based on a genetic algorithm to obtain the optimal parameter solution. Correspondingly, in the highly complex inversion problem of identifying groundwater pollution sources in three-dimensional fractured aquifers, although the preceding steps have constructed an inversion framework based on an integrated deep learning alternative model, transforming the physical problem into a mathematical optimization problem, the objective function typically exhibits high nonlinearity, non-convexity, and a large number of local extrema. Traditional gradient-based optimization methods are prone to getting trapped in local optima and cannot find a globally optimal parameter combination that truly reflects the groundwater pollution status within the vast parameter space. Furthermore, the search space is extremely complex because the parameters to be inverted include various heterogeneous variables such as pollution source location, release history, and hydrogeological parameters. Therefore, step S5 is introduced to perform global optimization based on a genetic algorithm to overcome the interference of local extrema and ultimately converge to the globally optimal parameter solution that minimizes the error between simulated and actual observations.

[0044] In one possible implementation, step S5 is performed as follows: First, constrained population initialization and initial fitness evaluation are conducted. The program first initializes an empty collection container, named Initial Population - Generation 0, to store the first generation of candidate solutions for the genetic algorithm. Then, the population size is set. The value is 30. The program enters an individual generation loop, which executes 30 times. In each iteration, a single individual is generated. .individual Essentially, it's a vector of decision variables, with dimensions corresponding to the number of parameters in the identified list of sensitive parameters, such as 14 dimensions. For an individual... Each gene in That is, decision variable components The program is within its corresponding prior knowledge boundary , Uniformly distributed random sampling is performed within the range of [450, 550] meters. For example, for the X coordinate of the pollution source, a value is randomly generated within the interval [450, 550] meters, such as 482.5; for matrix permeability, [the value is randomly generated within the interval of [450, 550] meters]. Sampling within the interval, such as This process ensures that all generated initial individuals are physically valid. The 30 generated individuals are then stored in the initial population set.

[0045] Next, the initial fitness batch evaluation is performed. The program iterates through every individual in the initial population. This calls the `fitness_function` encapsulated within the ensemble inversion framework. Inside the function, it first utilizes normalization parameters and the extreme feature values ​​saved during the training phase... The data is converted into a normalized vector in the [0,1] interval and then input into the trained alternative model. This alternative model uses its pre-trained weights and biases to predict the corresponding multi-well concentration time series within milliseconds. Subsequently, the root mean square error (RMSE) between the predicted sequence and the field-observed concentration data obtained in step S3 is calculated and converted into fitness values. The calculated 30 fitness values ​​are stored in an initial fitness score array. The program compares these scores, identifies the individual with the highest value, and assigns it to the globally optimal individual variable as a benchmark for subsequent evolution.

[0046] The algorithm then enters its core iterative evolutionary loop based on genetic operators. This loop will continue until a termination condition, such as the maximum number of iterations, is met. Or fitness convergence. In each generation cycle, the selection operation is performed first. Based on the fitness score of the current generation, a roulette wheel selection mechanism is used to select 30 individuals from the population to form a mating pool. Under this mechanism, each individual... The probability of being selected is directly proportional to its fitness. This means that individuals with high fitness (i.e., smaller error and closer to the actual pollution situation) have a greater chance of being selected and passing on their superior genes to the next generation, while individuals with low fitness are gradually eliminated.

[0047] Next, the crossover operation is performed. The program randomly selects pairs of parent individuals from the mating pool. Using crossover probability =0.9 determines whether to perform crossover. If the random number is less than 0.9, arithmetic crossover is performed on the decision variable vectors of this parent pair. Two offspring are generated using the formula: In the formula, It is a random variable uniformly distributed in the interval [0,1]. This formula shows that the offspring genes are a linear convex combination of the parent genes. For example, if the parents... The fracture density is 450, and the parent... The fracture density is 550, and =0.4, then the offspring The fracture density is 0.4 × 450 + 0.6 × 550 = 510. This operation allows the offspring to inherit the characteristic range of the parent and performs interpolation search in the solution space. The generated offspring are placed into a temporary next-generation population container.

[0048] The mutation operation is then performed, a crucial step in maintaining population diversity and preventing premature convergence. The program iterates through every individual in the next generation of the population. Every gene Based on the probability of mutation =0.2 determines whether mutation is performed. If mutation is performed, a non-uniform mutation strategy is used, and the calculation formula is as follows: In this formula, and This refers to the upper and lower bounds of the constraint set from which the gene originates; It is the index of the current evolutionary generation (e.g., generation 100). This is the preset total number of generations (500 generations). and It is a random number between [0,1]. The core of this formula lies in the adjustment factor. In the early stages of evolution The factor is relatively small, close to 1, and the variable-asynchronous length is relatively large, allowing the algorithm to perform a large-scale global exploration; as the evolution progresses... near As this factor approaches 0, the variable asynchronous length gradually decreases, and the algorithm switches to performing a fine-grained search within a local range. This dynamic adjustment strategy significantly improves the convergence accuracy of the algorithm.

[0049] After mutation, a new population is generated and evaluated. At this point, the individuals that have undergone crossover and mutation constitute the complete next-generation population. The program again calls the `fitness_function` to calculate the fitness of each new individual. Then, a global optimum update is performed: the individual with the highest fitness in the next-generation population is found and compared with the historical global optimum. If the individual with the highest fitness has a higher fitness, then that individual is used to update the global optimum. This ensures that the algorithm always retains the historical best solutions discovered during the search process.

[0050] Finally, convergence checks and optimal solution output are performed. At the end of each iteration, it is checked whether the maximum number of generations has been reached. The evolutionary cycle terminates if either condition is met, or if the variance of the fitness over several consecutive generations is less than a preset threshold. At this point, the decision variable vector stored in the globally optimal individual variable is the final optimal parameter solution. This solution contains a set of specific numerical values, such as the coordinates of the pollution source. =473.0, =102.5, =86.1, matrix permeability This set of optimal parameter solutions represents the combination of pollution source characteristics and geological parameters that best explains the current pollution phenomenon, given prior hydrogeological knowledge and field observation data.

[0051] In particular, during the inversion and optimization process for tracing groundwater pollution sources, traditional evaluation systems typically rely on the standard root mean square error (RMSE) as the sole criterion for evaluating the quality of candidate solutions. However, this single-dimensional evaluation method suffers from significant feature blindness when dealing with pollutant transport problems characterized by high nonlinearity and time-varying properties. Specifically, this method fails to deeply understand and utilize the unique physical relationships inherent in pollutant concentration time series data, namely the intrinsic connection between key event characteristics and global trends. The information value of pollutant concentration monitoring data is not uniformly distributed across all time points; the amplitude and timing of peak concentrations are crucial information characterizing the core physical features of pollution events, directly related to the maximum release intensity and peak period of the pollution source. The standard RMSE function treats errors at all time points equally, making it unable to distinguish between two solutions with fundamentally different qualities: one that predicts accurately near the peak point but performs poorly in other flat areas, and another that has a uniform overall error distribution but exhibits significant deviations at key peak points. This indiscriminate evaluation method results in the objective function lacking a deep understanding of the physical processes of pollution events. This not only leads to a mathematical loss of crucial information but also causes evaluation bias in a physical sense, failing to guide the optimization algorithm to converge to a comprehensive optimal solution. To address these technical shortcomings, this application introduces an adaptive fitness construction step based on key event features. The aim is to build a composite evaluation system that can intelligently identify and focus on key physical features, thereby guiding the global optimization algorithm's search process more accurately and efficiently.

[0052] In a possible preferred implementation, step S5, performing global optimization and parameter inversion on the integrated inversion framework based on a genetic algorithm to obtain the optimal parameter solution, includes: step S51, obtaining candidate solution vectors from the genetic algorithm; step S52, preprocessing the candidate solution vectors to obtain normalized candidate solution vectors; step S53, inputting the normalized candidate solution vectors into the trained alternative model to obtain the predicted concentration time series; step S54, calculating the fitness value based on the predicted concentration time series and the field observed concentration data; and step S55, performing global optimization and parameter inversion based on the fitness value to obtain the optimal parameter solution.

[0053] Step S51: During the iterative process of the genetic algorithm, each individual in the population represents a potential solution, and this individual is encoded as a vector containing multidimensional decision variables. The program extracts the gene sequence of the i-th individual from the current generation of the population, i.e., the candidate solution vector P. This vector covers all the key parameters to be inverted, including pollution source characteristic parameters and hydrogeological parameters. For example, if the number of decision variables is set to 14, the extracted candidate solution vector P may be [473.0, 102.5, 86.1, 0.95, 5.45 × 10]. -11[0.27, 430, 10.4, 10.5, 4.38, ...]. The first three values ​​represent the spatial coordinates of the pollution source (X, Y, Z, respectively, in meters); the fourth value represents the peak leakage rate of the sinusoidal release history (in mol / s); and the fifth value represents the matrix permeability (in meters). 2 The subsequent values ​​correspond to matrix porosity and fracture density (fractures / km) respectively. 3 The parameters include fracture dip angle (degrees), fracture orientation (degrees), and fracture aperture (mm). These values ​​are randomly generated within a preset prior knowledge boundary through crossover, mutation, and selection operations using a genetic algorithm.

[0054] Step S52: Due to the huge differences in the dimensions and orders of magnitude of the physical parameters in the candidate solution vector P, for example, the spatial coordinates are hundreds of meters, while the matrix permeability is only 10. -11 The magnitude of the input is so large that directly inputting it into the model can lead to instability in numerical calculations and weight bias. Therefore, it is necessary to utilize the statistical characteristics of the parameters (maximum values) determined during the training phase of the alternative model. and minimum value The program performs min-max normalization on vector P. It iterates through each component of vector P, linearly mapping it to the interval [0,1]. Taking the X coordinate of the pollution source as an example, if the prior search range is [450,550] and the current value is 473.0, then the normalized value is (473.0-450) / (550-450) = 0.23. Taking matrix permeability as an example, if the range is [4×10...]... -11 6×10 -11 The current value is 5.45 × 10 -11 The normalized value is (5.45-4) / (6-4) = 0.725. After this step, the original vector P with physical units is transformed into a dimensionless normalized candidate solution vector. Its format is adapted to the input layer requirements of deep learning models.

[0055] Step S53: The program calls the pre-trained deep learning alternative model with fixed parameters. The normalized vector obtained in step S52 is fed into the input feature tensor of the model. The model extracts local features through internal convolutional layers, captures long-term and short-term temporal dependencies through GRU and LSTM layers, and focuses on the weights of key time steps through an attention mechanism layer, performing forward propagation calculations. The output layer of the model generates the corresponding normalized concentration response, and restores it to a physically meaningful concentration value through an inverse normalization operation. The final output is the predicted concentration time series. This sequence contains concentration change data from all monitoring wells within a preset observation period. For example, for three monitoring wells (P1, P4, and P7) arranged parallel to the flow direction, the model outputs a 3×100 matrix (observation period of 1000 days, step size of 10 days). Specifically, the predicted sequence might show that well P1 reaches a peak concentration of 24.8 mol / m³ on day 400. 3 Well P7 reached its peak at 5.2 mol / m³ on day 800. 3 This results in a complete set of penetration curves that can reflect the geological and pollution source characteristics of the current candidate solutions, which can be used for subsequent comparison and evaluation with field observation data.

[0056] In a possible preferred embodiment, step S54, calculating the fitness value based on the predicted concentration time series and the field observed concentration data, includes: step S541, calculating the global consistency error between the predicted concentration time series and the field observed concentration data; step S542, quantifying the key event feature differences between the predicted concentration time series and the field observed concentration data to obtain a peak feature penalty term; and step S543, determining the fitness value based on the peak feature penalty term and the global consistency error.

[0057] Step S541 aims to quantify the degree of agreement between the two curves in terms of overall trend, ensuring that the inversion solution maintains physical rationality on a macroscopic level. The program calls the root mean square error formula for calculation, which is the same process as the objective function calculation in step S4 above, i.e. , Representing global consistency error, taking a real-world scenario with 3 monitoring wells (M=3) and 100 observation time points (T=100) as an example, if the model predicts that the concentration of monitoring well No. 1 at the 50th time step is 24.5 mol / m³, this would be considered a global consistency error. 3 The on-site observation value was 25.0 mol / m 3 The program will calculate the sum of squared residuals for this point and all other 299 data points, then take the average and square root. The calculated result... The value is 1.85, which represents the consistency error level of the current candidate solution in the global scope and serves as the basis for subsequent composite evaluation.

[0058] The specific implementation process of step S542 involves quantifying the key event feature differences between the predicted concentration time series and the field-observed concentration data to obtain a peak feature penalty term. This proactively identifies and quantifies peak features crucial for the qualitative characterization of pollution events by introducing domain knowledge. In a possible preferred implementation, step S542, quantifying the key event feature differences between the predicted concentration time series and the field-observed concentration data to obtain a peak feature penalty term, includes: quantifying the key event feature differences between the predicted concentration time series and the field-observed concentration data using the following formula: in, The peak feature penalty term, and These are weighted hyperparameters used to adjust the importance of peak amplitude and peak time error, respectively. and It is the maximum concentration value extracted from predicted concentration time series and field observation concentration data. and This is the index of the corresponding peak occurrence time. This refers to the time series step size, such as 10 days, used to standardize the units of measurement. The program first executes a peak detection algorithm on the input predicted and observed sequences to extract the maximum concentration values. and and its corresponding time index and Then, the extracted feature values ​​are substituted into the weighted penalty function. and These are preset weighting hyperparameters used to adjust the relative importance of amplitude error and time error. The settings of these two parameters are usually based on expert experience or the statistical characteristics of historical data. For example, if the study focuses more on the accuracy of source intensity, the weighting hyperparameters can be adjusted. Set to a larger value (e.g., 1.0); increase the value if you are more concerned with the accuracy of the leak occurrence time. Assume the peak value of the observed data is 25.0 mol / m. 3 It appeared on day 400 (index 40), while the predicted peak was 23.0 mol / m. 3 Appears on day 420 (index 42), set =1.0, =0.01, =10, then the peak feature penalty term =1.0×(23.0-25.0)^2+0.01×(10×(42-40))^2=8.0. This value intuitively quantifies the degree of failure of the candidate solution in reproducing the core plot of the pollution.

[0059] Step S543 aims to construct a composite evaluation index that balances globality and criticality. The program uses an adaptive fitness function to combine the results calculated in the first two steps. and Dynamic fusion is performed, and the formula is as follows: .in, It is a very small constant (e.g., 1×10⁻⁶). -6 This is to prevent the denominator from being zero. It is a balancing factor between 0 and 1, used to dynamically adjust the relative importance of global and local features. The setting of this value depends on the needs of the inversion stage: in the early stages of inversion, to ensure the algorithm quickly converges to a roughly correct solution space, a smaller value can be set. For example, a value of 0.3 emphasizes global fitting; in the later stages of inversion, to finely correct the peak shape, the value can be appropriately increased. Continuing with the example above, if... =1.85, =8.0, setting =0.4, then the denominator is calculated as (1-0.4)×1.85+0.4×sqrt{8.0}≈2.24, and the final fitness value is... ≈1 / 2.24≈0.446. The higher this fitness value, the more accurately the candidate solution captures the key pollution peak characteristics while ensuring the reasonableness of the global trend.

[0060] Step S55 involves global optimization and parameter inversion based on fitness values ​​to obtain the optimal parameter solution. This step is the core evolutionary stage of the genetic algorithm. The program evaluates and ranks all individuals in the population based on the fitness values ​​calculated in step S54. A roulette wheel selection mechanism is used, with individuals having a higher fitness value (e.g., 0.446) having a greater probability of being selected, thus passing on their superior genes (i.e., coordinates and hydrogeological parameters closer to the actual pollution source) to the next generation. Subsequently, a new offspring population is generated through single-point crossover (probability e.g., 0.9) and non-uniform mutation (probability e.g., 0.2). This process iterates until a termination condition is met (e.g., reaching the maximum number of iterations of 500 or the fitness no longer significantly increases). Finally, the decoded vector corresponding to the individual with the highest fitness is output, which is the optimal parameter solution. For example, the final output solution might include pollution source coordinates (X=514.3, Y=106.9, Z=84.8) and matrix permeability of 4.7 × 10⁻⁶. -11 m 2This solution not only minimizes the global error numerically, but also successfully reproduces the peak morphology and outbreak time of the pollution plume in a physical sense. The global optimization and parameter inversion process executed in step S55 includes roulette wheel selection based on fitness scores, crossover operation to generate offspring using the arithmetic crossover formula, mutation operation to maintain population diversity using a non-uniform mutation strategy, and the final logic for population update, convergence judgment, and optimal solution output. The specific implementation details and calculation formulas are completely consistent with the iterative evolution loop process based on genetic operators described in detail in the first embodiment of step S5 above, and will not be repeated here.

[0061] In step S6, the optimal parameter solution is comprehensively interpreted to obtain a pollution source characteristic report. This report includes the three-dimensional coordinates of the pollution source, time-series data representing the release history, and estimated values ​​of various hydrogeological parameters. In other words, although the previous steps have used a genetic algorithm to search for the globally optimal parameter solution in a complex parameter space, this solution is internally represented by a set of high-dimensional, abstract numerical vectors. For environmental engineers, decision-makers, or remediation teams, these cold numerical values ​​are difficult to directly translate into an intuitive understanding of the pollution site and effective engineering actions. These mathematically optimal solutions need to be restored to physically meaningful geological and pollution characteristics and presented in a visual manner to truly reveal the spatial location of the pollution source, leakage patterns, and key attributes of the aquifer. Therefore, step S6, which involves comprehensive interpretation of pollution source characteristics, transforms the abstract numerical results output by the algorithm into understandable and actionable engineering decision-making basis. By generating a comprehensive report containing three-dimensional coordinates, dynamic release history curves, and estimated hydrogeological parameters, it directly serves subsequent pollution risk assessment and remediation plan formulation.

[0062] In one possible implementation, step S6 is carried out as follows: First, parameter deconstruction is performed. The program reads the decision variable vector stored in the globally optimal individual. This vector is a one-dimensional array containing 14 components, each component corresponding to a specific physical parameter. Based on the parameter arrangement order defined in step S4, the program uses indexing operations to separate these values ​​one by one and assign them physical meanings. Based on the genetic algorithm inversion results, the deconstruction process is as follows: First, the first three components are extracted and identified as the three-dimensional spatial coordinates of the pollution source, for example, obtaining... =473.0 meters, =102.5 meters, =86.1 meters. This set of coordinates clarifies the specific location of the pollution source within the underground aquifer. Next, parameters characterizing the release history, namely the leakage rate feature value, are extracted. In this embodiment, such as the leakage rate... Parameterized as a sine function The deconstructed parameters determine the specific form of the function, thereby reconstructing the time-varying flux release sequence. Finally, estimated values ​​of key hydrogeological parameters, such as matrix permeability, are extracted. Crack density fracture dip angle and fracture aperture These deconstructed values ​​are stored as independent variables with definite physical units, preparing them for subsequent processing.

[0063] The results were then visualized. To make the deconstructed data more intuitive, a 3D visualization engine (such as Matplotlib, PyVista, or specialized groundwater simulation post-processing software) was used to convert the data into graphics. First, a 3D location map of the pollution source was drawn. In a 1000m×1000m×100m virtual model space, a transparent bounding box representing the study area was drawn, and a highlighted red sphere was marked at the corresponding position according to the deconstructed coordinates (473.0, 102.5, 86.1), visually showing that the pollution source is located on the left side of the model, near the bottom and slightly off-center. At the same time, the positions of the monitoring wells (P1, P4, P7) can be marked in the form of blue cylinders to form a visual comparison of spatial relative relationships. Next, a release history curve was drawn. The time t (0-1000 days) was used as the horizontal axis, and the deconstructed release flux was plotted as the vertical axis. Plot a continuous curve with (mol / s) as the vertical axis. For example... Figure 3 As shown, Figure 3 This is a release history curve from the deep learning-enhanced inversion-based groundwater pollution source localization method in fractured formations, according to an embodiment of this application. The curve exhibits a clear sine wave shape, with a peak value of approximately 1.0 mol / s and a period of approximately 1000 days. To demonstrate the reliability of the inversion, a 95% confidence interval shaded band (derived from genetic algorithm population statistics) can be superimposed on the graph. The narrower the shaded band, the lower the uncertainty of the inversion results. This visualization directly reveals whether the pollutant is continuously leaking or intermittently being released, and whether the current period is a peak leakage period.

[0064] Finally, the report is generated. The program automatically calls a document generation library (such as Python-docx or ReportLab) to integrate the quantified parameter values ​​and the generated images into a standardized pollution source characteristic report. The report content is structured: the first part is the pollution source location results, listing the inversion coordinates. The study includes a 3D location map with text indicating that the pollution source is located upstream of the study area; the second part is a reconstruction of the release history, showcasing... The graph and text description indicate that the pollution source experienced periodic fluctuations in release, reaching peak values ​​at specific time periods. The third part is the estimation of hydrogeological parameters, listing the comparison between the inverted values ​​and the prior range in tabular form, for example, indicating that the fracture density is 430 fractures / km². 3 Slightly below the prior average, this indicates a moderate degree of fracture development in the area, with permeability primarily controlled by the fracture aperture (4.38 mm). This report not only summarizes all technical details but also provides data-driven scientific interpretations, directly delivering them to the engineering team for subsequent borehole verification, groundwater remediation scheme design (such as well layout optimization), and long-term environmental risk management.

[0065] In summary, the deep learning-based enhanced inversion method for locating groundwater pollution sources in fractured strata, based on embodiments of this application, has been elucidated. Addressing the technical problems of high computational cost of numerical simulation, neglect of the heterogeneity of fractured aquifers, and low inversion accuracy due to insufficient utilization of time-series data in existing groundwater pollution source identification methods, this scheme first constructs a spatiotemporal concentration response dataset containing complex physical field characteristics based on hydrogeological and pollution source parameters. A deep learning alternative model is then trained based on this dataset to efficiently approximate the complex groundwater solute transport process, thus solving the problem of excessively long processing times in traditional numerical simulations. Furthermore, an integrated inversion framework is constructed, combining time-series concentration data from field observations with the alternative model, and using a genetic algorithm for global optimization. This not only handles nonlinear high-dimensional parameter spaces but also fully explores the spatiotemporal evolution patterns in the observation data, ultimately achieving simultaneous and accurate inversion of the three-dimensional coordinates of the pollution source, its dynamic release history, and key hydrogeological parameters. This effectively overcomes the shortcomings of traditional methods, such as low identification rate and poor timeliness under complex geological conditions.

[0066] Figure 4 This is a block diagram of a deep learning-based augmented inversion-based groundwater pollution source localization system for fractured strata, according to an embodiment of this application. Figure 4As shown, the deep learning-based enhanced inversion groundwater pollution source localization system 100 based on an embodiment of this application includes: a training module 110 and a localization module 120; the training module 110 includes: a spatiotemporal concentration response data construction unit 111, used to construct a spatiotemporal concentration response dataset based on a hydrogeological parameter set and a pollution source parameter set; and an alternative model training unit 112, used to train and validate an alternative model based on the spatiotemporal concentration response dataset to obtain a trained alternative model; the localization module 120 includes: an observation concentration data acquisition unit 121, used to acquire field observation concentration data; and an integrated... The inversion framework construction unit 122 is used to construct an integrated inversion framework based on field observation concentration data and a trained alternative model. The integrated inversion framework includes decision variables, objective function, and constraints. The optimal parameter solution analysis unit 123 is used to perform global optimization and parameter inversion on the integrated inversion framework based on a genetic algorithm to obtain the optimal parameter solution. The pollution source feature report generation unit 124 is used to perform comprehensive interpretation of pollution source features on the optimal parameter solution to obtain a pollution source feature report. The pollution source feature report includes the three-dimensional coordinates of the pollution source, time series data representing the release history, and estimated values ​​of various hydrogeological parameters.

[0067] Here, those skilled in the art will understand that the specific operations of each step in the above-described deep learning-based enhanced inversion system for locating groundwater pollution sources in fractured formations have been referenced above. Figures 1 to 3 The method for locating groundwater pollution sources in fractured strata based on deep learning-enhanced inversion has been described in detail, and therefore, its repeated description will be omitted.

Claims

1. A method for locating groundwater pollution sources in fractured strata based on deep learning-enhanced inversion, characterized in that, include: Training phase and positioning phase; The training phase includes: Based on the hydrogeological parameter set and the pollution source parameter set, a spatiotemporal concentration response dataset is constructed; The alternative model is trained and validated based on the spatiotemporal concentration response dataset to obtain the trained alternative model; The positioning phase includes: Obtain on-site concentration data; Based on field observation concentration data and trained alternative models, an integrated inversion framework is constructed, which includes decision variables, objective function and constraints. The optimal parameter solution is obtained by performing global optimization and parameter inversion based on genetic algorithm on the integrated inversion framework. The optimal parameter solution is subjected to comprehensive interpretation of pollution source characteristics to obtain a pollution source characteristic report, which includes the three-dimensional coordinates of the pollution source, time series data representing the release history, and estimated values ​​of various hydrogeological parameters.

2. The method for locating groundwater pollution sources in fractured strata based on deep learning-enhanced inversion according to claim 1, characterized in that, Based on hydrogeological parameter sets and pollution source parameter sets, a spatiotemporal concentration response dataset is constructed, including: A three-dimensional fracture network was randomly generated and the model was discretized from the set of hydrogeological parameters to obtain a discretized DFN / EPM model. Based on the pollution source parameter set, the physical field control equations of the discretized DFN / EPM model are solved and the global concentration field is generated to obtain the global spatiotemporal concentration field. Based on the virtual monitoring well location coordinate set, the spatiotemporal concentration field of the whole domain is used to extract the time series data of monitoring well concentration to obtain the multi-well concentration time series. Parameter-response pair correlations were performed on hydrogeological parameter sets, pollution source parameter sets, and multi-well concentration time series to obtain a spatiotemporal concentration response dataset.

3. The method for locating groundwater pollution sources in fractured strata based on deep learning-enhanced inversion according to claim 2, characterized in that, A three-dimensional fracture network was randomly generated from the hydrogeological parameter set, and the model was discretized to obtain a discretized DFN / EPM model, including: The discrete fracture network geometric entity set is obtained by randomly generating the discrete fracture network geometric entity set based on the set of hydrogeological parameters. A hybrid geometric model is constructed and a computational mesh is generated for the discrete fracture network geometric entity set to obtain a discretized DFN / EPM model.

4. The method for locating groundwater pollution sources in fractured strata based on deep learning-enhanced inversion according to claim 2, characterized in that, Based on the pollution source parameter set, the physical field control equations of the discretized DFN / EPM model are solved and the global concentration field is generated to obtain the global spatiotemporal concentration field, including: Based on the set of hydrogeological parameters, the discretized DFN / EPM model is coupled with the groundwater flow field solution to obtain the global water pressure distribution and the flow velocity field of the fracture domain and matrix domain. Time-domain pollutant migration is solved based on the velocity field of the fracture domain and matrix domain, hydrogeological parameter set and pollution source parameter set to obtain global concentration snapshot sequence; Data aggregation is performed on the global concentration snapshot sequence to obtain the global spatiotemporal concentration field.

5. The method for locating groundwater pollution sources in fractured strata based on deep learning-enhanced inversion according to claim 1, characterized in that, The alternative model includes a local temporal feature extraction module, a gated sequence encoding module, a long dependency modeling module, a key information weighting module, and a decoding module.

6. The method for locating groundwater pollution sources in fractured strata based on deep learning-enhanced inversion according to claim 1, characterized in that, To obtain the optimal parameter solution, a genetic algorithm-based global optimization and parameter inversion are performed on the ensemble inversion framework, including: Obtain candidate solution vectors from the genetic algorithm; The candidate solution vectors are preprocessed to obtain normalized candidate solution vectors; The normalized candidate solution vectors are input into the trained alternative model to obtain the predicted concentration time series; Fitness values ​​are calculated based on predicted concentration time series and field observation concentration data; Global optimization and parameter inversion are performed based on fitness values ​​to obtain the optimal parameter solution.

7. The method for locating groundwater pollution sources in fractured strata based on deep learning-enhanced inversion according to claim 6, characterized in that, Based on predicted concentration time series and field observation concentration data, fitness values ​​are calculated, including: Calculate the global consistency error between the predicted concentration time series and the field-observed concentration data; Key event feature differences were quantified between predicted concentration time series and field observation concentration data to obtain peak feature penalty terms; The fitness value is determined based on the peak feature penalty term and the global consistency error.

8. The method for locating groundwater pollution sources in fractured strata based on deep learning-enhanced inversion according to claim 7, characterized in that, To obtain a peak feature penalty term, the critical event feature difference is quantified between the predicted concentration time series and the field-observed concentration data. This includes: quantifying the critical event feature difference between the predicted concentration time series and the field-observed concentration data using the following formula: in, The peak feature penalty term, and These are weighted hyperparameters used to adjust the importance of peak amplitude and peak time error, respectively. and It is the maximum concentration value extracted from predicted concentration time series and field observation concentration data. and This is the index of the corresponding peak occurrence time. It is the step size of the time series.

Images