Underground water pollution risk prediction method fusing multi-source heterogeneous data

Abstract

The invention belongs to the technical field of groundwater environment, and relates to a multi-source heterogeneous data fused groundwater pollution risk prediction method, which comprises the following steps of: constructing a standardized multi-dimensional input data set containing a water quality state and an environment driving factor; constructing an initial hydrogeological topological graph of the underground water monitoring area; generating a dynamic physical state diagram reflecting the real-time hydraulic conduction capability; generating a node depth feature sequence capable of representing an independent change rule of each monitoring point; calculating a space source intensity attention matrix representing local emission abnormity; constructing a time conduction attention flow based on physical state constraint; and generating a global groundwater pollution risk prediction result, fusing the space source intensity attention matrix and the time conduction attention flow, and outputting the future risk probability of each node. According to the method, the problem that data adaptability and physical consistency of the model are difficult to consider due to the fact that a physical mechanism model and a data driving model are separated from each other in an existing prediction method is solved.

Description

Technical Field

[0001] This invention belongs to the technical field of groundwater environment and relates to a groundwater pollution risk prediction method that integrates multi-source heterogeneous data. Background Technology

[0002] Predicting groundwater pollution risks is crucial for environmental protection and water resource management; however, the prediction process faces numerous technical challenges. Groundwater systems are inherently heterogeneous, anisotropic, three-dimensional porous media, and their internal water flow and solute transport processes are influenced by complex hydrogeological conditions and dynamic external environmental factors. Furthermore, the monitoring data used to characterize this system typically exhibit spatial sparsity, multi-source heterogeneity, and temporal non-stationarity, making it difficult to construct reliable prediction models.

[0003] To address this issue, the industry typically employs two main technical approaches: one is based on physical process-based numerical simulation methods, such as the finite element method or finite difference method, which simulate the spatiotemporal distribution of pollutants by solving partial differential equations governing groundwater flow and solute transport. The other is a purely data-driven approach, utilizing time series analysis, statistical regression, or deep learning models to directly learn the intrinsic patterns of pollutant concentration changes from historical monitoring data. While numerical simulation methods can generally follow physical laws well, their accuracy heavily relies on the precise characterization of hydrogeological parameters. Obtaining these parameters is often costly and involves significant uncertainty, leading to limitations in parameter calibration and real-time data assimilation.

[0004] Based on the above problems, the present invention aims to solve the problem that the physical mechanism model and the data-driven model are separated in the existing prediction methods, making it difficult for the model to balance data adaptability and physical consistency. Summary of the Invention

[0005] To address the aforementioned problems, this invention provides a method for predicting groundwater pollution risk by integrating multi-source heterogeneous data.

[0006] A method for predicting groundwater pollution risk by fusing multi-source heterogeneous data includes the following steps:

[0007] S1. Acquire multi-source heterogeneous monitoring data, perform spatiotemporal alignment and cleaning on the multi-source heterogeneous monitoring data, and construct a standardized multidimensional input dataset containing water quality status and environmental driving factors.

[0008] S2. Construct an initial hydrogeological topology map of the groundwater monitoring area, abstract the monitoring wells as graph nodes, and initialize the edge weights between nodes according to static hydrogeological parameters.

[0009] S3. Generate a dynamic physical state diagram that reflects the real-time hydraulic conduction capacity, and dynamically correct the edge weights of the initial hydrogeological topology diagram using physical response rules.

[0010] S4. Extract the pollutant time-series characteristics of the monitoring well nodes and generate a node depth feature sequence that can characterize the independent change pattern of each monitoring point;

[0011] S5. Calculate the spatial source intensity attention matrix that characterizes local emission anomalies and identify monitoring nodes that exhibit high-risk source characteristics at the current moment;

[0012] S6. Construct a time-propagating attention flow based on physical state constraints, and use physical mechanisms to guide the propagation path of data features in the graph network;

[0013] S7. Generate prediction results for groundwater pollution risks across the entire region, integrate spatial source strength attention matrix and temporal transmission attention flow, and output the future risk probability of each node.

[0014] A further aspect of this invention involves constructing a standardized multidimensional input dataset that includes water quality status and environmental driving factors, comprising the following steps:

[0015] Acquire time-series data on groundwater pollutant concentrations collected by the monitoring well network, as well as environmental driving data including rainfall intensity and surface evaporation collected by meteorological monitoring stations;

[0016] Acquire remote sensing image data covering the monitoring area, analyze the surface soil type and vegetation cover, and use them as static environmental parameters;

[0017] Pollutant concentration time series data are mapped to a preset time window, and environmental driving data and static environmental parameters are spatially indexed and aligned with the geographic coordinates of monitoring wells to generate a standardized multidimensional input dataset.

[0018] A further aspect of this invention involves constructing an initial hydrogeological topographic map of a groundwater monitoring area, comprising the following steps:

[0019] Define monitoring wells as graph nodes, calculate the Euclidean distance between monitoring wells, and establish a connection relationship between nodes if the Euclidean distance is less than a preset spatial connection threshold.

[0020] Obtain the static permeability coefficient of the aquifer, and calculate the initial hydraulic connection strength of each connection relationship based on the ratio of the static permeability coefficient to the Euclidean distance.

[0021] Assigning the initial hydraulic connection strength as the initial weight of the graph edges generates an initial hydrogeological topology map that serves as a static skeleton of physical constraints.

[0022] A further aspect of the present invention generates a dynamic physical state diagram reflecting real-time hydraulic conduction capability, comprising the following steps:

[0023] Call the pre-set parameterized physical response function library, which contains parameterized formulas describing the disturbance of the groundwater flow field by surface infiltration;

[0024] Input the environmental driving data and static environmental parameters into the parameterized physical response function library to calculate the hydraulic conduction enhancement coefficient matrix under the current time window;

[0025] The edge weights of the initial hydrogeological topology map are modulated in real time using the hydraulic conduction enhancement coefficient matrix to generate a dynamic physical state map.

[0026] A further aspect of this invention utilizes a hydraulic conduction enhancement coefficient matrix to modulate the edge weights of an initial hydrogeological topology map in real time, comprising the following steps:

[0027] Extract rainfall intensity from environmental driving data and soil type from static environmental parameters;

[0028] The infiltration potential coefficient is determined based on soil type, and the hydraulic conduction enhancement coefficient is calculated using a parameterized formula in conjunction with rainfall intensity.

[0029] The initial weights of the initial hydrogeological topology map are multiplied element by element with the corresponding hydraulic conduction enhancement coefficients to obtain dynamic edge weights, thus forming a dynamic physical state map.

[0030] A further aspect of this invention involves generating a node depth feature sequence capable of characterizing the independent variation patterns of each monitoring point, comprising the following steps:

[0031] Extract time-series data of groundwater pollutant concentrations from a standardized multidimensional input dataset and input them into a temporal convolutional network or a long short-term memory network.

[0032] Temporal convolutional networks or long short-term memory networks are used to extract the concentration fluctuation trend, periodic changes and abrupt change signals of each monitoring well in the time dimension;

[0033] The extracted signals are mapped into high-dimensional vectors to generate node depth feature sequences corresponding to each graph node.

[0034] A further aspect of this invention involves calculating the spatial source intensity attention matrix characterizing local emission anomalies, comprising the following steps:

[0035] Perform a linear transformation on the node depth feature sequence to generate query vectors and key vectors;

[0036] Calculate the dot product of the query vector and the key vector to obtain the feature similarity between nodes;

[0037] The feature similarity is normalized using Softmax to generate a spatial source-strength attention matrix that focuses on high-concentration emission sites.

[0038] A further aspect of this invention involves constructing a time-propagating attention flow based on physical state constraints, comprising the following steps:

[0039] By referencing the dynamic physics state graph, an edge weight matrix containing real-time physical meaning is extracted;

[0040] The edge weight matrix is ​​used as a structured bias term for the attention mechanism and superimposed on the original attention score calculated based on the node depth feature sequence;

[0041] The superimposed scores are normalized to generate a time-transmitted attention stream that characterizes the probability of pollutant migration with the groundwater flow field.

[0042] A further aspect of this invention uses the edge weight matrix as a structured bias term in the attention mechanism, comprising the following steps:

[0043] Introduce learnable balance parameters to the edge weight matrix;

[0044] The weighted edge weight matrix is ​​added to the original attention score generated by the node depth feature sequence, thus forcing the features to propagate along a physically feasible hydraulic transmission path.

[0045] A further aspect of this invention generates a comprehensive groundwater pollution risk prediction result, comprising the following steps:

[0046] A gated fusion unit is used to weight and fuse the spatial source-strength attention matrix and the temporal transmission attention stream to generate physical-data dual-constraint features;

[0047] The physical-data dual-constraint features are input into the fully connected prediction layer to decode and obtain the predicted pollutant concentration values ​​of each monitoring well node and the blank area of ​​the monitoring well. Spatial interpolation is then performed based on the physical connection relationship to obtain the predicted pollutant concentration values ​​of the blank area of ​​the monitoring well.

[0048] By comparing the predicted pollutant concentration with the preset risk threshold, a groundwater pollution risk prediction result is generated, which includes the location of the pollution source and the range of migration and diffusion.

[0049] In summary, the present invention has the following beneficial technical effects:

[0050] 1. By constructing a dynamic physical state diagram to constrain the time-transmitted attention flow, the physical impact of real-time environmental driving factors on the groundwater flow field can be quantified and integrated into the feature propagation process of the data-driven model. This dynamic physical state diagram adjusts the edge weights representing the hydraulic connection strength between monitoring well nodes in real time based on parameters such as rainfall intensity and soil type. This weight matrix then serves as a structured bias term for the attention mechanism, guiding feature information to aggregate only along physically reasonable paths. This mechanism avoids the model learning spatial correlations that are physically impossible due to data randomness.

[0051] 2. By parallel computing of the spatial source intensity attention matrix and the temporal transmission attention stream, and finally fusing them using a gating mechanism, the system can distinguish and integrate the local generation characteristics and spatial transmission characteristics of pollutants. The spatial source intensity attention is calculated based on the abnormal patterns of the node's own temporal data, without relying on physical connections, and can effectively identify newly emerging and independent pollution emission sources. The temporal transmission attention stream, under the constraints of the physical map, simulates the migration process of pollutants from upstream to downstream. The final gating fusion unit dynamically determines the weights of the two features based on the data. This mechanism enables the system to perform a deeper analysis of the causes of pollution events, realize the location of pollution sources, and track the migration path of pollution plumes.

[0052] 3. Utilizing temporal convolutional networks to extract the temporal features of pollutants from monitoring well nodes can capture the multi-scale dynamic patterns of data over time. Temporal convolutional networks, through stacked causal-perforated convolutional layers, exponentially expand the receptive field, extracting long-term trends, periodic fluctuations, and short-term abrupt changes in the original concentration sequence in parallel across different time scales, and encoding these dynamic patterns into high-dimensional node depth feature sequences. Compared to traditional time series analysis methods, this mechanism can more fully exploit the complex nonlinear relationships inherent in time series data, providing more representative input features for subsequent attention calculations and risk prediction, and improving the model's ability to characterize the dynamic evolution of different types of pollution events. Attached Figure Description

[0053] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the accompanying drawings used in the description of the embodiments or the prior art will be briefly introduced below. The drawings are used to provide a further understanding of the present invention.

[0054] Figure 1 A flowchart illustrating an embodiment of this application is disclosed.

[0055] Figure 2 Structural schematic diagrams of embodiments of this application are disclosed. Detailed Implementation

[0056] The following is in conjunction with the appendix Figure 1 - Figure 2 A preferred description of the present invention is provided below.

[0057] See attached document Figure 1 This invention proposes a groundwater pollution risk prediction method that integrates multi-source heterogeneous data, comprising the following steps:

[0058] S1. Acquire multi-source heterogeneous monitoring data, perform spatiotemporal alignment and cleaning on the multi-source heterogeneous monitoring data, and construct a standardized multidimensional input dataset containing water quality status and environmental driving factors.

[0059] S2. Construct an initial hydrogeological topology map of the groundwater monitoring area, abstract the monitoring wells as graph nodes, and initialize the edge weights between nodes according to static hydrogeological parameters.

[0060] S3. Generate a dynamic physical state diagram that reflects the real-time hydraulic conduction capacity, and dynamically correct the edge weights of the initial hydrogeological topology diagram using physical response rules.

[0061] S4. Extract the pollutant time-series characteristics of the monitoring well nodes and generate a node depth feature sequence that can characterize the independent change pattern of each monitoring point;

[0062] S5. Calculate the spatial source intensity attention matrix that characterizes local emission anomalies and identify monitoring nodes that exhibit high-risk source characteristics at the current moment;

[0063] S6. Construct a time-propagating attention flow based on physical state constraints, and use physical mechanisms to guide the propagation path of data features in the graph network;

[0064] S7. Generate prediction results for groundwater pollution risks across the entire region, integrate spatial source strength attention matrix and temporal transmission attention flow, and output the future risk probability of each node.

[0065] In one embodiment of the present invention, step S1 includes the following steps:

[0066] Acquire time-series data of groundwater pollutant concentrations collected by the monitoring well network, as well as corresponding environmental driving data collected by meteorological monitoring stations. The environmental driving data includes at least rainfall intensity and surface evaporation. Acquire remote sensing image data covering the monitoring area, and analyze the surface soil type and vegetation cover as static environmental parameters.

[0067] Pollutant concentration time series data are mapped to a preset time window, and environmental driving data and static environmental parameters are spatially indexed and aligned with the geographic coordinates of monitoring wells to generate a standardized multidimensional input dataset.

[0068] Specifically, the data processing terminal performs the following steps. It asynchronously acquires raw data from multiple distributed data sources and generates structured input for subsequent models through a series of standardized preprocessing steps. The terminal periodically polls the sensor data acquisition systems deployed in the monitoring well network via database interface protocols such as ODBC or JDBC to obtain raw pollutant concentration time-series data streams containing unique well identifiers, high-frequency timestamps, and concentrations of various pollutants such as nitrates, nitrites, and heavy metal ions. The terminal also calls meteorological data services via HTTP API to obtain raw environmental-driven data streams from various meteorological monitoring stations associated with the monitoring area. These streams include timestamps, rainfall intensity, and surface evaporation. Finally, the terminal downloads in batches a collection of satellite remote sensing imagery covering the entire monitoring area, such as multispectral imagery from Sentinel-2 or the Landsat series, via FTP or an interface accessing cloud storage object services.

[0069] The data processing terminal performs spatiotemporal alignment processing on the acquired multi-source heterogeneous monitoring data. In the time dimension, the system presets a unified time window T, such as 1 hour, and employs a time resampling algorithm to aggregate the original pollutant concentration time-series data stream and the original environmental driving data stream into this time window. Typically, the mean value within the window is used to smooth high-frequency noise and achieve uniform temporal resolution. In the spatial dimension, the data processing terminal processes the remote sensing image dataset, using reflectance data from the near-infrared and red bands to calculate the Normalized Difference Vegetation Index (NDVI) to quantify vegetation cover. Multispectral data is then input into a pre-trained supervised classification model, such as a random forest model, to identify the surface soil type and assign a category code to each pixel, thereby generating a rasterized static environmental parameter layer.

[0070] Using the geographic coordinates of each monitoring well as an index, the data processing terminal precisely extracts the vegetation cover value and soil type code for each monitoring well location from the static environmental parameter layer. For environmentally driven data where the location does not completely overlap with the monitoring well, spatial interpolation algorithms, such as the inverse distance weighting (IDW) method, are used to estimate the equivalent rainfall intensity and surface evaporation at each monitoring well location based on the distance and observation values ​​between the monitoring well and its surrounding meteorological stations. After data alignment, data cleaning and standardization are performed. The data processing terminal fills in missing values ​​in the aligned data using linear interpolation or spline interpolation algorithms and employs a method based on... Criteria or local anomaly factor algorithms identify and eliminate potential measurement outliers. To eliminate the influence of different dimensions between features, a min-max normalization method is used to linearly map all numerical features to the interval between 0 and 1.

[0071] The water quality status data, i.e. pollutant concentration, is combined with environmental driving factor data, including rainfall intensity, surface evaporation, soil type coding, and vegetation cover, at each time window to construct a logical three-dimensional tensor. The dimensions are the number of monitoring wells, time step, and feature dimension, respectively. This tensor is the standardized multidimensional input dataset, which serves as the unified input for the subsequent fusion of physical mechanisms and data-driven models.

[0072] Vegetation cover is quantified by calculating the Normalized Difference Vegetation Index (NDVI), and its calculation formula is expressed as follows:

[0073]

[0074] NIR represents the reflectance value in the near-infrared band, and R represents the reflectance value in the red band. Both are dimensionless parameters, and their values ​​are determined by the quantization level of the remote sensing image sensor.

[0075] Among them, multi-source heterogeneous monitoring data refers to various data types with different sources, structures, and spatiotemporal resolutions, including point-like time-series data collected from monitoring wells, point-like time-series data collected from meteorological stations, and areal spatial data collected from remote sensing satellites. Environmentally driven data refers to external environmental factors that can affect groundwater processes and pollutant migration, mainly rainfall intensity and surface evaporation, usually measured in mm / h and mm / day.

[0076] Supervised classification models are trained on labeled sample data and can map new input data to predefined categories, such as mapping the spectral features of remote sensing images to soil types like sand and clay. Static environmental parameters refer to surface features that remain stable or change slowly over a longer timescale, mainly surface soil type and vegetation cover. The former is usually represented by dimensionless category codes, while the latter are continuous values ​​between 0 and 1.

[0077] The time window is the baseline duration for time alignment. The time window is set according to the needs of the application scenario. For example, if the application scenario is a chemical industrial park that requires a rapid response to sudden leakage events, the time window T can be set to 1 hour; while in the long-term trend analysis of regional seasonal changes, it can be set to 24 hours.

[0078] For example, suppose the system needs to process data from monitoring well W01 between 10:00 and 11:00. The data processing terminal obtains two pollutant concentration records for W01 during this period: 0.5 mg / L at 10:15 and 0.7 mg / L at 10:45. The system also obtains a rainfall intensity of 5 mm / h recorded at 10:30 from a nearby weather station. The set time window T is 1 hour. During the time alignment phase, the system performs mean aggregation on the concentration data of W01 from 10:00 to 11:00, obtaining an average concentration of (0.5 + 0.7) / 2 = 0.6 mg / L for this time window. The rainfall intensity data is directly included in this time window. The terminal retrieves the latest remote sensing image covering the geographic coordinates of W01, resolving the near-infrared reflectance (NIR) of the pixel at that location to be 0.5 and the red reflectance (R) to be 0.1. Based on this, the system calculates the vegetation coverage characterization value NDVI = (0.5-0.1) / (0.5+0.1)≈0.67.

[0079] The pre-trained soil classification model identifies the location as sandy loam and assigns a category code of 3. These features are combined into an original feature vector, such as [0.6, 5, 0.67, 3]. In the standardization stage, assuming that based on historical data, pollutant concentrations range from 0 to 2 mg / L and rainfall intensity ranges from 0 to 20 mm / h, the vector is subjected to min-max normalization. The concentration term becomes (0.6-0) / (2-0) = 0.3, and the rainfall intensity term becomes (5-0) / (20-0) = 0.25. Vegetation cover itself is between 0 and 1 or close to this range and can be left unprocessed or processed according to the actual range. Soil type codes, as category features, are usually one-hot encoded or left unchanged. After complete processing, a standardized feature vector corresponding to monitoring well W01 in the 10:00-11:00 time window is generated in the multidimensional input dataset, for example, [0.3, 0.25, 0.67, 3], and stored in the corresponding position of the three-dimensional tensor.

[0080] In one embodiment of the present invention, step S2 includes the following steps:

[0081] The monitoring wells are defined as nodes in the graph. Based on the Euclidean distance between the monitoring wells and the static permeability coefficient of the aquifer, the connection relationship between the nodes is established. The initial hydraulic connection strength of each connection relationship is calculated and assigned as the initial weight of the graph edge to generate an initial hydrogeological topology map, which serves as the static skeleton of the physical constraints.

[0082] Specifically, this process is executed by the graph construction module. The graph construction module aims to transform discrete monitoring well sites into a physically meaningful topological structure. It traverses the monitoring well list processed in step S1, creating a graph node object in memory for each well's unique identifier. The module retrieves the 3D geographic coordinates (x, y, z) of each well from the associated geographic information database and stores them as core attributes in the corresponding graph node object. The module establishes connections between nodes to form graph edges. It employs a fully connected traversal strategy, calculating the Euclidean distance between any two monitoring well nodes i and j. This calculation result is compared to a preset spatial connectivity threshold. If the calculated Euclidean distance is less than or equal to the threshold, it is determined that the two monitoring wells have a potential direct hydrogeological connection, and an undirected edge is established between the corresponding graph nodes. Conversely, if the Euclidean distance is greater than the threshold, it is considered that the spatial distance between the two is too great, the hydraulic connection is weak, or there is no direct connection, and no edge is established.

[0083] For each established edge, the graph construction module calculates its initial weight, defined as the initial hydraulic connection strength. This calculation follows a simplified version of Darcy's Law, which states that the hydraulic connection strength is directly proportional to the permeability of the aquifer and inversely proportional to the length of the flow path. The graph construction module retrieves the static permeability coefficient values ​​of the aquifers corresponding to the regions containing nodes i and j from a static hydrogeological parameter database. This static permeability coefficient value is divided by the Euclidean distance between nodes i and j, and the result is assigned as the initial weight of the edge connecting these two nodes. After calculating the weights of all edges, the system integrates all graph nodes, edges, and their corresponding initial weights to form a weighted undirected graph data structure. This data structure is stored in the form of an adjacency matrix or adjacency list, thus constituting the initial hydrogeological topology graph. This graph provides the static physical framework for subsequent dynamic adjustments and feature propagation.

[0084] Euclidean distance between any two monitoring well nodes i and j The calculation is as follows:

[0085]

[0086] Initial weights of the edge connecting node i and node j The initial hydraulic connection strength is calculated using the following formula:

[0087]

[0088] Among them, Euclidean distance It is a scalar that measures the straight-line distance between two nodes in three-dimensional space, where... and These are the three-dimensional geographic coordinates of nodes i and j, respectively. The static permeability coefficient of the aquifer. Physical parameters characterize the ability of porous underground media to allow fluid to pass through, typically measured in m / d. This value is considered static because it depends on geological structures and does not change in the short term. Its values ​​are derived from geological exploration reports or obtained through inversion from pumping test data. Graph nodes are digital abstractions of physical monitoring wells and are the basic units of the graph.

[0089] During the calculation, Take the arithmetic or harmonic mean of the permeability coefficients at corresponding locations of node i and node j. Initial hydraulic bond strength. It is an indicator that quantifies the groundwater exchange capacity between two monitoring wells under undisturbed conditions, and its dimensions are: In a physical sense, it can be understood as a measure of the fluid exchange rate under a unit head difference.

[0090] The initial hydrogeological topology map is a mathematical graph composed of a set of nodes, edges, and weights, describing the basic hydraulic conduction network structure within the monitoring area under static conditions. The spatial connectivity threshold is set based on the geological scale of the monitoring area and the average connectivity range of the hydrogeological units. Assuming the application scenario is an alluvial fan aquifer in a plain area with good connectivity, the threshold can be set to 1000m to 3000m.

[0091] For example, suppose the system needs to construct an initial hydrogeological topology map for three monitoring wells W01, W02, and W03. The graph construction module obtains the geographic coordinates, assuming the coordinates are on the same plane, namely W01 (100, 200), W02 (400, 600), and W03 (1200, 200), in meters. The system sets the spatial connectivity threshold to 1000 meters, and the graph construction module calculates the Euclidean distance between each pair of nodes. The distance between W01 and W02 is calculated. for m. Distance between W01 and W03 for m. Distance between W02 and W03 for m.

[0092] The calculated distance is compared with the threshold, because m and m, the system establishes edges between W01 and W02, and between W02 and W03, and because No edge is established between W01 and W03. The system retrieves the static permeability coefficient of the aquifer in this area from the hydrogeological parameter database. The weight is 15m / d. The system calculates the initial weights of the edges W01 and W02. The weights between W02 and W03 This generates an initial hydrogeological topology map, which contains three nodes W01, W02, and W03, and two edges (W01, W02) and (W02, W03), with weights of 0.03 and 0.0168, respectively.

[0093] In one embodiment of the present invention, step S3 includes the following steps:

[0094] Call the preset parametric physical response function library, which contains parametric formulas describing the disturbance law of surface infiltration on the groundwater flow field; input the environmental driving data and static environmental parameters of step S1 into the parametric physical response function library to calculate the hydraulic conduction enhancement coefficient matrix under the current time window.

[0095] The edge weights of the initial hydrogeological topology map in step S2 are modulated in real time using the hydraulic conduction enhancement coefficient matrix to generate a dynamic physical state map that instantly reflects changes in the groundwater flow field and the migration paths of potential pollutants.

[0096] Specifically, the dynamic graph generation module executes within each time window to generate a dynamic physical state graph reflecting the real-time hydraulic conduction capacity under the current environmental conditions. This module loads the initial hydrogeological topology map generated in step S2 from the system memory and obtains its static edge weight matrix. The dynamic graph generation module extracts the environmental driving data corresponding to the current time window from the standardized multidimensional input dataset generated in step S1, specifically the rainfall intensity at each monitoring well location and static environmental parameters, i.e., the surface soil type.

[0097] The dynamic graph generation module calls a pre-built parametric physical response function library. This library stores multiple software collections describing response functions for different physical processes. The module selects a parametric formula from this library to characterize the disturbance of groundwater flow field by surface infiltration. The module takes extracted rainfall intensity and surface soil type data as input and feeds them into the selected parametric formula. For each edge in the initial hydrogeological topology map, the module combines the environmental parameters corresponding to the two nodes connected by the edge and calculates the hydraulic conduction enhancement coefficient value using the formula. This coefficient value quantifies the enhancement factor of groundwater conduction capacity along the path under the current rainfall infiltration conditions. Combining the calculation results of all edges generates a hydraulic conduction enhancement coefficient matrix corresponding to the graph structure.

[0098] The module performs a weight modulation operation, multiplying the static edge weight matrix of the initial hydrogeological topology map element-wise with the newly generated hydraulic conduction enhancement coefficient matrix. The result of this operation is a new dynamic edge weight matrix, in which the weight of each edge reflects the physical state changes caused by the current environmental driving factors in real time. The module uses this new dynamic edge weight matrix to replace the original static weights, thereby generating a dynamic physical state map that reflects the changes in the groundwater flow field and the migration paths of potential pollutants in real time. This map, as the output, provides a physically reasonable constraint for the subsequent temporal characteristics of pollutant propagation.

[0099] For the edge connecting node i and node j, its dynamic weight The calculation method is as follows:

[0100]

[0101] Hydraulic conduction enhancement coefficient The result is calculated using a parameterized physical response function, with an exemplary function form as follows:

[0102]

[0103] The hydraulic conduction enhancement coefficient matrix is ​​a matrix of the same dimension as the graph adjacency matrix, with elements... It is a dimensionless scalar with a value greater than or equal to 1, representing the degree of enhancement of hydraulic connection. It is the updated dynamic edge weight, and its dimensions are the same as the initial weight. Same, for . These are the initial weights derived from the initial hydrogeological topology map. It is the effective rainfall intensity affecting edge (i, j), which can be the average of the rainfall intensity corresponding to nodes i and j, in mm / h. It is the soil type that affects edge (i, j), and can be a comprehensive representative type of the soil at two nodes.

[0104] It is an infiltration potential coefficient related to soil type, and its value is obtained from a preset lookup table based on soil type. The search returned results, such as those corresponding to clay. The values ​​are relatively low, ranging from 0.1 to 0.3, while those for sandy soil are lower. The values ​​are relatively high, ranging from 0.7 to 0.9. It is an empirical saturation rate constant, with units of 1000 m / s. This value is used to control the rate at which the enhancement effect tends to saturate as rainfall intensity increases. It is calibrated based on historical data and is typically in the range of 0.05-0.2.

[0105] The Parametric Physical Response Function Library is a code library containing multiple mathematical models, each describing a specific physical process through adjustable parameters. For example, in addition to surface infiltration models, it may also include models describing the effects of river replenishment or artificial water injection.

[0106] For example, the initial hydrogeological topology map generated in step S2 contains edges (W01, W02) and (W02, W03), with initial weights of respectively. and Assuming the current time window, the system obtains the following rainfall intensity from the standardized multidimensional input dataset: 5 mm / h at location W01, with a soil type of sandy loam (code 3); 8 mm / h at location W02, with a soil type of sandy loam (code 3); and 8 mm / h at location W03, with a soil type of clay (code 4).

[0107] The infiltration potential coefficient of sandy loam soil is found in the system's preset lookup table. It is 0.8, the clay It is 0.2, and the saturation rate constant is 0.2. Set as Calculate the hydraulic conduction enhancement coefficient for edge (W01, W02). The soil connected to this side is all sandy loam, therefore Effective rainfall intensity mm / h. Calculate the coefficients of edge (W02, W03). The edge connects sandy loam and clay, and their average infiltration potential coefficient is taken. Effective rainfall intensity mm / h. Therefore, .

[0108] Calculate the dynamic weights. . Output a dynamic physical state graph with its topology unchanged, but the weights of edges (W01, W02) and (W02, W03) have been updated to 0.041472, which better reflects the current impact of rainfall infiltration. and 0.0214284 .

[0109] In one embodiment of the present invention, step S4 includes the following steps:

[0110] Input the time-series data of groundwater pollutant concentration from step S1 into a temporal convolutional network or a long short-term memory network; extract the concentration fluctuation trend, periodic changes, and abrupt change signals of each monitoring well in the time dimension; map the extracted signals into high-dimensional vectors to generate a node depth feature sequence corresponding to each graph node.

[0111] Specifically, this is performed by the temporal feature extraction module. This module aims to extract high-dimensional abstract features that characterize the inherent dynamic laws from the original pollutant concentration time series. From the standardized multidimensional input dataset generated in step S1, this module extracts the corresponding pollutant concentration time series data for each monitoring well node. The system slices the time series data of each node using a sliding window, forming multiple input sequence segments of length L, where L is the preset historical time step. The temporal feature extraction module feeds these input sequence segments in batches into a pre-trained temporal convolutional network. This temporal convolutional network consists of multiple stacked residual blocks, each containing two causal dilated convolutional layers, a weight normalization layer, and an activation function layer.

[0112] Causal convolution ensures that the output at time step t depends only on the input at and before time t, conforming to the inherent logic of time series. Dilated convolution, through an exponentially increasing dilation factor, allows the network's receptive field to expand rapidly with increasing layer depth, effectively capturing long-term dependencies in time series, such as slow concentration accumulation trends or seasonal cyclical fluctuations, without adding excessive computational burden. When an input sequence of length L passes through this network, each convolutional kernel scans the sequence at different time scales, extracting hierarchical features layer by layer, from local minor fluctuations to global trend changes.

[0113] The module takes the output of the last residual block of the temporal convolutional network at the last time step. This output is a fixed-dimensional vector, which is a highly condensed and abstract representation of the entire input sequence segment. The numerical combinations within it already contain complex information such as concentration fluctuation trends, periodic changes, and abrupt signals from the original sequence. For each monitoring well node, the system generates such a high-dimensional vector through this process. The set of these vectors constitutes the node depth feature sequence corresponding to each graph node.

[0114] The core operation of temporal convolutional networks is causal dilated convolution. For a one-dimensional input sequence s and a convolutional kernel F, the convolution operation at time step t can be represented as:

[0115]

[0116] in, This represents historical data points in the input sequence that are sampled at intervals relative to the current time point t, based on the inflation factor d. It represents the weight of the convolution kernel at position i; It is the size of the convolution kernel; It is the expansion factor, which is usually set to a power of 2 at different levels of the network, such as 1, 2, 4, 8, to achieve rapid expansion of the receptive field. The formula clearly reflects causality because when calculating the output at time t, the summation only involves the input points at time t and before.

[0117] Temporal convolutional networks (TCNNs) are deep learning models specifically designed for processing sequential data. They capture temporal dependencies through hierarchical convolutional structures and offer better parallel computing capabilities and more stable gradients compared to recurrent neural networks (RNNs). Long Short-Term Memory (LSTM) networks are another alternative sequence model. By introducing input, forget, and output gates, they control the flow and retention of information and are also adept at handling long-term dependencies. Node deep feature sequences are high-dimensional real-valued vectors, representing the result of encoding the original temporal data through a deep neural network, mapping the dynamic patterns of the time series onto the feature space.

[0118] For example, suppose the system processes pollutant concentration data from monitoring well W01, setting the input sequence length L to 12 time steps, and the generated node depth feature sequence dimension to 8. The temporal feature extraction module obtains the standardized concentration values ​​of W01 for the most recent 12 time windows from the standardized multidimensional input dataset, forming the input sequence, for example [0.30, 0.31, 0.32, 0.33, 0.35, 0.36, 0.55, 0.40, 0.38, 0.37, 0.36, 0.35]. This sequence exhibits a pattern of first slowly rising, then a sudden spike, and then falling back. This input sequence is fed into a pre-trained temporal convolutional network. The shallow convolutional layers in the temporal convolutional network may mainly identify small changes between adjacent time points in the sequence, while the deep, perforated convolutional layers, due to their large receptive field, can simultaneously observe the growth trend at the beginning of the sequence and the abrupt events in the middle.

[0119] After computation and nonlinear transformations across all layers, the temporal convolutional network outputs an 8-dimensional vector at the 12th time step. Assuming this output vector is [0.41, -0.25, 0.88, 0.12, -0.56, 0.33, 0.05, -0.19], this vector represents the node depth feature sequence of W01 at that time step. A high value in one dimension of this vector, such as 0.88 in the third dimension, may represent a strong response of the network to a recent abrupt signal event. This process is applied in parallel to all other monitoring well nodes, generating a unique depth feature sequence for each node.

[0120] In one embodiment of the present invention, step S5 includes the following steps:

[0121] Calculate the spatial source intensity attention matrix to characterize local emission anomalies and identify monitoring nodes that exhibit high-risk source characteristics at the current moment.

[0122] Self-attention calculation is performed on the node depth feature sequence generated in step S4 to generate query vectors and key vectors; the feature similarity between nodes is calculated, and abnormally high value nodes are identified based solely on the intensity of the data features themselves, without relying on graph structure information; through Softmax normalization, a spatial source intensity attention matrix focusing on high concentration emission points is generated.

[0123] Specifically, the spatial attention analysis module performs this task. Its core mission is to identify monitoring nodes that exhibit significant pollution source characteristics within the current time window using a data-driven approach. The input to this module is the set of node depth feature sequences generated for all monitoring wells in step S4. This module employs a self-attention mechanism, passing the node depth feature sequence of each monitoring well node in parallel through three independent linear transformation layers with different learnable parameters to generate a query vector Q, a key vector K, and a value vector V. These three vectors may be the same as or different from the input feature sequence in terms of dimension and play different roles. The query vector represents the exploration signal actively emitted by the node, and the key vector represents the feature identifier of the node itself that can be matched.

[0124] The spatial attention analysis module calculates the feature similarity between nodes. For any node's query vector Q, the module calculates its dot product with the key vectors K of all nodes, including itself. This dot product operation produces a scalar value, which quantifies the similarity between two nodes in the deep feature space. The results of the dot products of all query vectors and all key vectors are summarized to form the original attention score matrix. To ensure numerical stability and prevent gradient vanishing or exploding during model training, the module divides each element in this matrix by a scaling factor, typically set to the square root of the key vector dimension. To transform the original attention scores into a probabilistic weight distribution, the module independently applies the Softmax function to normalize each row of the scaled attention score matrix.

[0125] After processing with the Softmax function, the element values ​​in each row are between 0 and 1, and the sum of all elements in the row is 1. The resulting matrix is ​​the spatial source intensity attention matrix. In this matrix, if the values ​​in a certain column are generally high, it means that the monitoring well node corresponding to that column possesses depth characteristics. For example, a drastic concentration change makes it the focus of attention for most other nodes at the feature level. Such high-attention nodes, identified solely by the strength of the data features without relying on any pre-set graph structure information, are identified as high-risk anomaly sources at the current moment. Calculation of the spatial source intensity attention matrix. The core formula is as follows:

[0126]

[0127] The query vector Q, key vector K, and value vector V are obtained by multiplying the input node depth feature sequence matrix X by three trainable weight matrices. , and What was obtained, that is , , . and The dimension is Where N is the number of monitoring well nodes, It is the dimension of the key vector. It involves performing matrix multiplication of the query matrix Q and the transpose of the key matrix K to obtain... The original attention score matrix. It is a scaling factor used to adjust the magnitude of the score. The setting is usually based on model design and experimental results, with typical values ​​of 64 or 128. Spatial source strength attention matrix yes A matrix whose elements This represents the attention weight of node i to node j, which reflects the degree of influence of node j as a potential source of pollution on node i.

[0128] Self-attention computation is a neural network mechanism that allows the model to dynamically evaluate the mutual importance between different parts of a dataset as it processes the input data. The Softmax function acts on each row of the matrix, converting the scores into a probability distribution.

[0129] For example, suppose we are currently processing three monitoring wells, W01, W02, and W03. As known from step S4, W01 exhibits an abnormal pattern in its node depth feature sequence due to a sudden concentration change. Assume that after linear transformation, the resulting two-dimensional bond vector... =[1.5, 2.0] (representing a strong feature of W01). =[0.2, 0.3] (representing the stable characteristics of W02). =[0.3, 0.2] (representing the stable features of W03). Query vector for W01. =[1.4, 1.8]. Key vector dimension The scaling factor is Calculate the raw attention score of W01 for all nodes.

[0130] and The dot product is . and The dot product is . and The dot product is Thus, the original score row vector corresponding to W01 is [5.7, 0.82, 0.78]. Scaling using a scaling factor yields [5.7 / 1.414, 0.82 / 1.414, 0.78 / 1.414]. [4.03, 0.58, 0.55]. Apply the Softmax function to this row vector to calculate... , , The values ​​are approximately 56.26, 1.78, and 1.73, respectively.

[0131] The normalized weights are 56.26 / (56.26+1.78+1.73). 0.941, 1.78 / 59.77 0.029, 1.73 / 59.77 0.029. In the first row of the generated spatial source strong attention matrix, the resulting vector is approximately [0.941, 0.029, 0.029]. The results show that W01's attention is mostly focused on itself, and in the attention calculation of other nodes, the column corresponding to W01 also receives a high score due to its unique key vector, thus highlighting its characteristics as a high-risk source in the entire matrix.

[0132] In one embodiment of the present invention, step S6 includes the following steps:

[0133] Using the dynamic physical state graph generated in step S3, extract the edge weight matrix containing real-time physical meaning; use the edge weight matrix as a structured mask or bias term of the attention mechanism to force the aggregation process of the node depth feature sequence in step S4, allowing features to propagate only along physically feasible hydraulic conduction paths; calculate the correlation strength between nodes after physical constraints, and generate a time-conducting attention flow characterizing the probability of pollutants migrating with the groundwater flow field.

[0134] Specifically, this is executed by the Physically Constrained Attention module, whose purpose is to construct a data feature propagation path that enforces physical laws. The Physically Constrained Attention module receives two sets of core inputs: one is the dynamic physical state graph generated in step S3 within the current time window, from which the module directly extracts the edge weight matrix containing real-time physical meaning; the other is the set of node depth feature sequences generated for all monitoring wells in step S4. This module inputs the node depth feature sequence set into an independent linear transformation layer to generate a new query vector Q' and key vector K' for each node. The Physically Constrained Attention module calculates the pure data-driven original attention score matrix. This calculation process is similar to step S5, that is, by calculating the dot product of the query vector Q' of each node and the key vector K' of all nodes, and then scaling it, a score matrix that initially reflects the feature similarity between nodes is obtained.

[0135] Unlike step S5, this module introduces physical constraints at this crucial step. The module uses the edge weight matrix extracted from the dynamic physics graph as a structured bias term, directly adding it to the previously calculated data-driven original attention score matrix. In this operation, if two nodes are physically connected in the dynamic physics graph, their corresponding positive weight values ​​will increase the attention score between these two nodes; conversely, if there is no physical connection between the two nodes (i.e., the weight is zero), their attention score will not receive any physical enhancement. This design forcibly injects prior knowledge of the physical mechanism into the formation process of attention weights.

[0136] The module normalizes the hybrid score matrix, modulated by the physical bias term, using the Softmax function. The Softmax function converts the modulated scores into a probability distribution, generating the final attention weight matrix. This matrix is ​​the time-transferred attention stream, where each weight value simultaneously embodies both the similarity at the data feature level and the feasibility of the physical hydraulic transport path, characterizing the probability of pollutants migrating from one node to another under current hydrogeological conditions. (Time-transferred attention stream) The calculation formula is as follows:

[0137]

[0138] in, and The query and key matrix is ​​obtained by performing a new linear transformation on the node depth feature sequence. It is a real-time edge weight matrix extracted from the dynamic physical state graph, and its elements are... This is dynamic weight. . It is a learnable scalar parameter used to balance the relative importance of data-driven attention scores and physical constraint bias terms. This parameter is automatically optimized during model training using the backpropagation algorithm.

[0139] Time-based attention flow yes A matrix whose elements This represents the degree of attention node i pays to the features of node j under physical constraints, which is also the relative probability of pollutants being transmitted from node j to node i.

[0140] Structured masks or bias terms are techniques that incorporate graph structure information into attention mechanisms. Using a mask usually means setting the score of unconnected locations to negative infinity so that their probability is zero, while using a bias term means adding the physical connection strength as an additional scoring term.

[0141] For example, the physics-constrained attention module receives a dynamic physics state graph with its edge weight matrix. It is known in the middle , ,and The module simultaneously receives the node depth feature sequences of W01, W02, and W03, generating new queries and key vectors. Suppose we now need to calculate the attention weights of node W02 regarding other nodes, i.e. The second row of the matrix. Assume the query vector for W02. The value is [0.8, 0.9], while the bond vectors of W01, W02, and W03 are... The values ​​are [0.9, 0.8], [1.0, 1.1], and [-0.5, 0.1], respectively, representing the key vector dimensions. The scaling factor is 1.414. Calculate the raw, purely data-driven score: W02's score against W01 is... W02's score for himself was W02's score against W03 was... The original score row vector is [1.018, 1.266, -0.220].

[0142] Introducing physical biases, assuming the trained parameters The physical bias term is Adding the two together yields the modulated score: [1.018+0.415, 1.266+0, -0.220+0.214] = [1.433, 1.266, -0.006]. Applying the Softmax function to [1.433, 1.266, -0.006], the final attention weight row vector is obtained, approximately [0.449, 0.382, 0.169]. This result represents the portion of the temporally transmitted attention stream corresponding to the W02 node, indicating that W02's attention is primarily allocated to itself and the upstream node W01, which has a direct strong physical connection.

[0143] In one embodiment of the present invention, step S7 includes the following steps:

[0144] By employing a cross-attention mechanism or gated fusion unit, the spatial source intensity attention matrix and the temporal transmission attention flow are weighted and fused to generate physical-data dual-constraint features. The physical-data dual-constraint features are then input into a fully connected prediction layer to decode and obtain the predicted pollutant concentration values ​​for each monitoring well node and the blank area of ​​the monitoring well. The predicted values ​​are then compared with a preset risk threshold to generate groundwater pollution risk prediction results that include the location of the pollution source and the range of migration and diffusion.

[0145] Specifically, this is executed by the fusion prediction module, which is responsible for integrating the attention information of different dimensions generated in the preceding steps and decoding and outputting the final risk prediction result. This module receives the spatial source-strength attention matrix from step S5 and the temporal transmission attention stream from step S6. This module applies a linear transformation layer independent of the generated query and key vectors to convert the node deep feature sequence generated in step S4 into a value vector V. The module utilizes a gated fusion unit to dynamically integrate the source-strength features and the physical transmission features; this unit performs the two information aggregation operations in parallel.

[0146] The first operation multiplies the spatial source-strength attention matrix with the value vector V matrix to obtain the source-strength aggregation feature matrix. Each row of this matrix represents the global anomaly source information aggregated by a node from the data level. The second operation multiplies the temporal transmission attention stream with the same value vector V matrix to obtain the transmission aggregation feature matrix. Each row of this matrix represents the upstream transmission information aggregated by a node along a physically feasible path. For each node, the gated fusion unit concatenates the two aggregation feature vectors and inputs the concatenated long vector into a fully connected layer containing a Sigmoid activation function to generate a gated scalar between 0 and 1. This gated scalar serves as a dynamic weight, used to perform a weighted summation of the source-strength aggregation feature and the transmission aggregation feature to generate a physical-data dual-constraint feature.

[0147] The module inputs the physical-data dual-constraint features generated for each node into a fully connected prediction layer for decoding. This prediction layer typically consists of one or two neural network layers with non-linear activation functions such as ReLU. Its task is to map the high-dimensional fused feature vector into scalar values, i.e., the predicted pollutant concentration values ​​for the corresponding monitoring well node in one or more future time steps. This module not only generates predicted values ​​for nodes with existing monitoring wells, but also performs interpolation calculations on virtual nodes in the monitoring area that are not equipped with monitoring wells, based on their physical connection relationships with surrounding nodes, to achieve concentration prediction with full coverage.

[0148] The system compares the predicted pollutant concentrations from all nodes with preset risk thresholds one by one. If a node's predicted value exceeds the threshold, the node is marked as high-risk. The system then aggregates all high-risk nodes and their predicted concentration values, combining this with geographic location information to generate a visualized groundwater pollution risk prediction result. This result is typically presented as a risk level map or warning list, identifying potential pollution source locations and the migration and diffusion range of the pollution plume over a future period. The gating fusion unit generates physical-data dual-constraint features. The calculation process is as follows:

[0149] First, calculate the source strength aggregation feature. and transmission aggregation features :

[0150]

[0151]

[0152] Then, calculate the gating weight G:

[0153]

[0154] Finally, the fused features are generated:

[0155]

[0156] Pollutant concentration prediction values Calculated from the fully connected prediction layer:

[0157]

[0158] Among them, the physical-data dual constraint feature It is a node feature representation that integrates data-driven source orientation information and physical mechanism-driven path transmission information. It is a spatial source strong attention matrix. It is time that transmits the flow of attention. Through linear transformation The resulting value matrix, where X is the node depth feature sequence matrix. This indicates a vector concatenation operation. It is the Sigmoid activation function. and These are the learnable weights and biases of the gating layer. This indicates element-wise multiplication. The function mapping represents the fully connected prediction layer.

[0159] The risk threshold is set based on national or local groundwater quality standards and anomaly judgment criteria based on the statistical distribution of historical monitoring data, such as using the 95th percentile of the normal background concentration as the threshold. The gating fusion unit is a neural network structure that learns dynamic gating signals to determine how to combine information streams from different sources.

[0160] For example, suppose the system needs to predict the future contamination risk of node W02, and the module receives the spatial source strong attention matrix. and the flow of attention over time Suppose the value vectors V of the three nodes are respectively , , .

[0161] Calculate the source strength aggregation characteristics of W02 ,need The second row, assuming it's [0.4, 0.5, 0.1], then... Calculate the conduction-aggregation characteristics of W02. Using the data obtained in S6 The second line [0.449, 0.382, 0.169], then Calculate the gate value. Concatenate the two features to obtain... Assuming the gate layer parameters and The calculated result is then passed through the Sigmoid function to obtain the gate value. .

[0162] Computational Physics-Data Dual Constraints The fused feature is input into the prediction layer, which is assumed to map it to a standardized concentration prediction value of 0.65. If the set risk threshold is 0.5, since 0.65 > 0.5, the system determines that node W02 has a high pollution risk and includes it in the groundwater pollution risk prediction results, marking it as an early warning state.

[0163] See appendix Figure 2 This invention also proposes a groundwater pollution risk prediction system that integrates multi-source heterogeneous data, comprising the following modules:

[0164] The multi-source data standardization module is used to acquire multi-source heterogeneous monitoring data, perform spatiotemporal alignment and cleaning on the multi-source heterogeneous monitoring data, and construct a standardized multidimensional input dataset that includes water quality status and environmental driving factors.

[0165] The static hydrogeological topology module is used to construct the initial hydrogeological topology map of the groundwater monitoring area, abstract the monitoring wells as graph nodes, and initialize the edge weights between nodes according to static hydrogeological parameters.

[0166] The Dynamic Physical State Diagram module is used to generate a dynamic physical state diagram that reflects the real-time hydraulic conduction capacity, and dynamically corrects the edge weights of the initial hydrogeological topology map using physical response rules.

[0167] The node feature extraction module is used to extract the temporal features of pollutants from the monitoring well nodes and generate a node depth feature sequence that can characterize the independent change patterns of each monitoring point.

[0168] The source strength attention module is used to calculate the spatial source strength attention matrix that characterizes local emission anomalies and identify monitoring nodes that exhibit high-risk source characteristics at the current moment.

[0169] The physical constraint propagation module is used to construct a time-propagated attention flow based on physical state constraints, and to guide the propagation path of data features in the graph network using physical mechanisms;

[0170] The risk fusion and prediction module is used to generate groundwater pollution risk prediction results across the entire region. It integrates the spatial source strength attention matrix and the temporal transmission attention flow to output the future risk probability of each node.

[0171] Each of the modules can be implemented in whole or in part through software, hardware, or a combination thereof. It supports hardware embedded in or independent of the processor in the computer device, and also supports software stored in the memory of the computer device, so that the processor can call and execute the operations corresponding to each of the above modules.

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

Claims

1. A method for predicting groundwater pollution risk by integrating multi-source heterogeneous data, characterized in that, Includes the following steps: S1. Acquire multi-source heterogeneous monitoring data, perform spatiotemporal alignment and cleaning on the multi-source heterogeneous monitoring data, and construct a standardized multidimensional input dataset containing water quality status and environmental driving factors. S2. Construct an initial hydrogeological topology map of the groundwater monitoring area, abstract the monitoring wells as graph nodes, and initialize the edge weights between nodes according to static hydrogeological parameters. S3. Generate a dynamic physical state diagram that reflects the real-time hydraulic conduction capacity, and dynamically correct the edge weights of the initial hydrogeological topology diagram using physical response rules. S4. Extract the pollutant time-series characteristics of the monitoring well nodes and generate a node depth feature sequence that can characterize the independent change pattern of each monitoring point; S5. Calculate the spatial source intensity attention matrix that characterizes local emission anomalies and identify monitoring nodes that exhibit high-risk source characteristics at the current moment; S6. Construct a time-propagating attention flow based on physical state constraints, and use physical mechanisms to guide the propagation path of data features in the graph network; S7. Generate prediction results for groundwater pollution risks across the entire region, integrate spatial source strength attention matrix and temporal transmission attention flow, and output the future risk probability of each node.

2. The groundwater pollution risk prediction method based on multi-source heterogeneous data according to claim 1, characterized in that, Constructing a standardized multidimensional input dataset that includes water quality status and environmental driving factors involves the following steps: Acquire time-series data on groundwater pollutant concentrations collected by the monitoring well network, as well as environmental driving data including rainfall intensity and surface evaporation collected by meteorological monitoring stations; Acquire remote sensing image data covering the monitoring area, analyze the surface soil type and vegetation cover, and use them as static environmental parameters; Pollutant concentration time series data are mapped to a preset time window, and environmental driving data and static environmental parameters are spatially indexed and aligned with the geographic coordinates of monitoring wells to generate a standardized multidimensional input dataset.

3. The groundwater pollution risk prediction method based on multi-source heterogeneous data according to claim 1, characterized in that, Constructing an initial hydrogeological topographic map of the groundwater monitoring area includes the following steps: Define monitoring wells as graph nodes, calculate the Euclidean distance between monitoring wells, and establish a connection relationship between nodes if the Euclidean distance is less than a preset spatial connection threshold. Obtain the static permeability coefficient of the aquifer, and calculate the initial hydraulic connection strength of each connection relationship based on the ratio of the static permeability coefficient to the Euclidean distance. Assigning the initial hydraulic connection strength as the initial weight of the graph edges generates an initial hydrogeological topology map that serves as a static skeleton of physical constraints.

4. The groundwater pollution risk prediction method based on multi-source heterogeneous data according to claim 1, characterized in that, Generating a dynamic physical state diagram reflecting real-time hydraulic conduction capacity includes the following steps: Call the pre-set parameterized physical response function library, which contains parameterized formulas describing the disturbance of the groundwater flow field by surface infiltration; Input the environmental driving data and static environmental parameters into the parameterized physical response function library to calculate the hydraulic conduction enhancement coefficient matrix under the current time window; The edge weights of the initial hydrogeological topology map are modulated in real time using the hydraulic conduction enhancement coefficient matrix to generate a dynamic physical state map.

5. The groundwater pollution risk prediction method based on multi-source heterogeneous data according to claim 4, characterized in that, The edge weights of the initial hydrogeological topology map are modulated in real time using the hydraulic conduction enhancement coefficient matrix, including the following steps: Extract rainfall intensity from environmental driving data and soil type from static environmental parameters; The infiltration potential coefficient is determined based on soil type, and the hydraulic conduction enhancement coefficient is calculated using a parameterized formula in conjunction with rainfall intensity. The initial weights of the initial hydrogeological topology map are multiplied element by element with the corresponding hydraulic conduction enhancement coefficients to obtain dynamic edge weights, thus forming a dynamic physical state map.

6. The groundwater pollution risk prediction method based on multi-source heterogeneous data according to claim 1, characterized in that, Generating node depth feature sequences that can characterize the independent variation patterns of each monitoring point includes the following steps: Extract time-series data of groundwater pollutant concentrations from a standardized multidimensional input dataset and input them into a temporal convolutional network or a long short-term memory network. Temporal convolutional networks or long short-term memory networks are used to extract the concentration fluctuation trend, periodic changes and abrupt change signals of each monitoring well in the time dimension; The extracted signals are mapped into high-dimensional vectors to generate node depth feature sequences corresponding to each graph node.

7. The groundwater pollution risk prediction method based on multi-source heterogeneous data according to claim 1, characterized in that, Calculating the spatial source-intensity attention matrix characterizing local emission anomalies includes the following steps: Perform a linear transformation on the node depth feature sequence to generate query vectors and key vectors; Calculate the dot product of the query vector and the key vector to obtain the feature similarity between nodes; The feature similarity is normalized using Softmax to generate a spatial source-strength attention matrix that focuses on high-concentration emission sites.

8. The groundwater pollution risk prediction method based on multi-source heterogeneous data according to claim 1, characterized in that, Constructing a time-propagating attention flow based on physical state constraints includes the following steps: By referencing the dynamic physics state graph, an edge weight matrix containing real-time physical meaning is extracted; The edge weight matrix is ​​used as a structured bias term for the attention mechanism and superimposed on the original attention score calculated based on the node depth feature sequence; The superimposed scores are normalized to generate a time-transmitted attention stream that characterizes the probability of pollutant migration with the groundwater flow field.

9. The groundwater pollution risk prediction method based on multi-source heterogeneous data according to claim 8, characterized in that, Using the edge weight matrix as a structured bias term in the attention mechanism includes the following steps: Introduce learnable balance parameters to the edge weight matrix; The weighted edge weight matrix is ​​added to the original attention score generated by the node depth feature sequence, thus forcing the features to propagate along a physically feasible hydraulic transmission path.

10. The groundwater pollution risk prediction method fusion method according to claim 1, characterized in that, Generating a comprehensive groundwater pollution risk prediction result includes the following steps: A gated fusion unit is used to weight and fuse the spatial source-strength attention matrix and the temporal transmission attention stream to generate physical-data dual-constraint features; The physical-data dual-constraint features are input into the fully connected prediction layer to decode and obtain the predicted pollutant concentration values ​​of each monitoring well node and the blank area of ​​the monitoring well. Spatial interpolation is then performed based on the physical connection relationship to obtain the predicted pollutant concentration values ​​of the blank area of ​​the monitoring well. By comparing the predicted pollutant concentration with the preset risk threshold, a groundwater pollution risk prediction result is generated, which includes the location of the pollution source and the range of migration and diffusion.

Images