Knowledge graph-based pollution site tracing method and decision system

Abstract

This invention relates to the field of artificial intelligence technology, specifically to a knowledge graph-based method and decision-making system for tracing contaminated sites. The method includes the following steps: analyzing the flow direction deviation angle and water flow velocity to generate effective components; fusing these components with topographic gradients to generate a dynamic coupling coefficient; using water level characteristics to compensate for this coefficient to obtain relational edge weights and construct migration relational edges; establishing a knowledge graph through network topology splicing; performing optimization comparison to extract candidate paths; fusing spatial attenuation features to obtain an index to evaluate path confidence and generate a decision. In this invention, the distribution network of underground pollution migration resistance is finely reconstructed through multi-source dynamic feature coupling, deeply integrating the macroscopic graph structure with the microscopic hydrological dynamic attenuation mechanism. This completely overcomes the limitations of single historical comparison and surface interpolation inference, achieving reliable optimization and ranking of topological nodes for all source tracing paths, and effectively improving the objectivity and authenticity of the environmental traceability decision-making system.

Description

Technical Field

[0001] This invention relates to the field of artificial intelligence technology, and in particular to a method and decision-making system for tracing the source of contaminated sites based on knowledge graphs. Background Technology

[0002] Artificial intelligence (AI) is a comprehensive technology system centered on computer systems simulating human intelligent activities. It encompasses core aspects such as data representation and storage, knowledge modeling and reasoning, machine learning and pattern recognition, natural language understanding and generation, and decision optimization. By performing structured representation and correlation analysis on multi-source data, AI enables the cognition, analysis, and decision support for complex problems. It has been widely applied in scenarios such as environmental monitoring, industrial control, and smart cities. Knowledge graphs, as an important branch of AI, construct entity relationship networks to semantically represent discrete data, providing fundamental support for causal analysis and source tracing in complex systems.

[0003] Traditional knowledge graph-based methods and decision-making systems for tracing contaminated sites involve constructing entity nodes using pollutant monitoring data, geospatial information, and historical emission records. A graph structure is formed by defining the relationships between pollution sources, migration paths, and receptors. In practice, spatial interpolation is used to generate distribution maps from collected soil sample concentrations. Suspected pollution sources are screened based on pollutant types and concentration thresholds. Migration path connections are established by combining groundwater flow data and terrain slope information. Simultaneously, rule-based matching is used to correlate and compare historical enterprise emission records with current monitoring results. Pollution sources are inferred based on pre-defined causal chains. During the decision-making process, nodes are back-analyzed according to the chronological order of pollutant diffusion paths, and comparisons are made with similar pollution patterns in an existing case library to form a source tracing judgment.

[0004] Traditional methods rely solely on spatial interpolation and rule matching to construct spectral structures and on surface concentration thresholds to screen suspected pollution sources. In the face of complex hydrogeological environments, they cannot accurately depict the dynamic characteristics and migration impedance of groundwater. The inference mechanism based on chronological backtracking and comparison of similar cases lacks quantitative assessment support of the physical decay law of pollutants. This results in structural biases in the connection of source tracing paths, making it difficult to guarantee the objective facts of source tracing in complex and compound pollution scenarios. This leads to serious consequences such as unclear pollution liability and errors in environmental governance decisions. Summary of the Invention

[0005] To address the technical problems of traditional methods that rely solely on spatial interpolation and rule matching to construct graph structures, and on surface concentration thresholds to screen suspected pollution sources, which fail to accurately depict groundwater dynamics and migration impedance in complex hydrogeological environments, and whose inference mechanisms based on chronological backtracking and comparison of similar cases lack quantitative assessment support of the physical decay laws of pollutants, resulting in structural biases in the connection of source tracing paths, making it difficult to guarantee the objective facts of source tracing determination in complex and complex pollution scenarios, and leading to serious consequences such as unclear pollution liability and erroneous environmental governance decisions, this invention provides a pollution site source tracing method based on knowledge graphs.

[0006] To achieve the above objectives, this invention employs a knowledge graph-based method for tracing the source of contaminated sites, comprising the following steps: S1: Collect the slope angle, groundwater flow direction angle and water flow velocity of the contaminated site. Based on the slope angle, perform spatial direction difference extraction on the groundwater flow direction angle to obtain the flow direction deviation angle. Based on the flow direction deviation angle, perform direction vector projection analysis on the water flow velocity to generate effective water flow components. S2: Call the effective water flow component to obtain the slope and spatial distance of adjacent land units, extract the spatial change rate of the slope based on the spatial distance, obtain the terrain gradient, perform dynamic feature coupling processing on the terrain gradient and the effective water flow component, and generate a dynamic coupling coefficient. S3: Call the dynamic coupling coefficient, collect the benchmark water level and real-time groundwater level of the contaminated site, perform water level compensation adjustment on the dynamic coupling coefficient, obtain the relation edge weights, perform feature mapping between adjacent land parcel units based on the relation edge weights, and construct pollution migration relation edges; S4: Perform network topology splicing processing on the pollution migration relationship edge and adjacent land parcel units to construct a contaminated site knowledge graph. Input the contaminated site knowledge graph into the Dijkstra algorithm for optimization and comparison to generate candidate source tracing paths. S5: Perform spatial decay feature fusion on the path node spacing of the candidate source tracing path and the preset pollutant decay rate to obtain the path decay index. Based on the path decay index, perform confidence assessment and ranking screening on the candidate source tracing path to generate a source tracing decision scheme.

[0007] As a further aspect of the present invention, the effective water flow component includes the actual migration driving component after coupling the groundwater flow direction with the slope direction, the dynamic coupling coefficient includes hydraulic gradient, permeability coefficient and hydraulic conductivity, the relation edge weight is specifically the pollution migration association intensity value between plots after water level compensation correction, the contaminated site knowledge graph includes plot nodes, migration relation edges and edge weight distribution, and the source tracing decision scheme includes candidate path priority ranking, path confidence score and pollution source pointing result.

[0008] As a further aspect of the present invention, the specific steps of S1 are as follows: S101: Collect the horizontal elevation slope angle, groundwater flow direction angle and instantaneous water flow velocity of the test points of the contaminated site. Extract the projection component of the horizontal elevation slope angle by spatial orthogonal decomposition. Perform Hadamard product operation on the groundwater flow direction angle and the projection component to obtain the site angle distribution matrix. S102: Perform discrete partial derivative calculation on adjacent node elements in the site angle distribution matrix to extract feature difference sequence, sort the elements of the feature difference sequence by absolute value comparison, remove abnormal discrete values ​​at both ends with a preset ratio, and perform inverse trigonometric transformation of vector angle on the sorted elements to generate flow direction deviation angle; S103: Substitute the flow direction deviation angle into the cosine function to extract the cosine mapping coefficient value, perform scalar dot product calculation on the instantaneous water flow velocity and the cosine mapping coefficient value to extract the projected velocity scalar set, and perform sliding average smoothing conversion on the projected velocity scalar set according to the preset sliding window length to generate effective water flow components.

[0009] As a further aspect of the present invention, the specific steps of S2 are as follows: S201: Collect the topographic slope and lateral spatial distance of adjacent plots in the contaminated site, perform discrete difference on the topographic slope to extract the node slope difference value sequence, divide the node slope difference value sequence by the lateral spatial distance to extract spatial variation features, and obtain the spatial variation rate matrix. S202: Perform grid mapping on the spatial rate of change matrix to extract the spatial gradient direction vector, obtain the elevation parameter of the contaminated site, perform scalar product calculation on the spatial gradient direction vector and the elevation parameter to extract the initial gradient scalar set, sort the elements in the initial gradient scalar set, and establish the terrain gradient. S203: Based on the effective water flow component and the terrain gradient, perform tensor product operation to extract the dynamic distribution feature matrix, obtain the fluid density constant of the contaminated site, multiply the dynamic distribution feature matrix with the fluid density constant to extract the absolute kinetic energy index set and perform normalization scaling to generate the dynamic coupling coefficient.

[0010] As a further aspect of the present invention, the step of performing discrete difference extraction of node slope difference value sequence for terrain slope refers to taking the terrain slope of adjacent plots in the contaminated site as input, calculating the slope difference value at each node in the spatial grid by subtracting the slope of adjacent nodes, and performing the calculation of the slope difference value sequentially along the extension direction of the lateral spatial distance. The arrangement order of each element in the node slope difference value sequence corresponds one-to-one with the arrangement order of nodes of adjacent plots in the lateral space.

[0011] As a further aspect of the present invention, the specific steps of S3 are as follows: S301: Collect the benchmark water level and real-time groundwater level of the contaminated site, perform subtraction calculation based on the benchmark water level and real-time groundwater level to extract the drop sequence, perform division calculation based on the drop sequence and benchmark water level to extract the fluctuation ratio set, perform exponential transformation on the fluctuation ratio set, and establish the water level compensation coefficient. S302: Multiply the dynamic coupling coefficient and the water level compensation coefficient to extract the connectivity dynamic parameter, perform range standardization mapping on the connectivity dynamic parameter, extract the interval feature value, obtain the connection damping constant, and perform difference operation based on the interval feature value and the connection damping constant to obtain the relation edge weight; S303: Call the relation edge weights, extract the spatial topological coordinate matrix of adjacent plot units, concatenate the relation edge weights with the topological coordinate matrix to extract the associated feature vectors, sort the associated feature vectors in descending order and remove the end nodes, assign the remaining connections with corresponding parameters, and generate pollution migration relation edges.

[0012] As a further aspect of the present invention, the specific steps of S4 are as follows: S401: Call the pollution migration relationship edge, collect spatial topological feature parameters, perform connection reconstruction based on the spatial topological feature parameters to extract the topological connection matrix, define node entities according to the topological connection matrix, aggregate the node entities with the pollution migration relationship edge, and establish a knowledge graph of the contaminated site. S402: Call the contaminated site knowledge graph, read historical test data to extract contaminated nodes and target nodes, retrieve the connection paths between nodes according to the topology connection matrix, accumulate the contaminated migration relationship edge association parameters step by step to extract the cumulative impedance value, and establish a path impedance weight set. S403: Call the path impedance weight set, sort the values ​​in ascending order and extract the first value, define the path corresponding to the first value as the benchmark connected path, compare the remaining values ​​with the preset impedance threshold to remove items exceeding the threshold, obtain the set of qualified connected paths, and concatenate the paths in the set of qualified connected paths with the benchmark connected paths to extract the optimized related subset and generate candidate traceability paths.

[0013] As a further aspect of the present invention, the preset impedance threshold is obtained by extracting the cumulative impedance values ​​corresponding to all known pollution migration paths in historical test data, performing a difference operation between the maximum and minimum values ​​among all cumulative impedance values, multiplying the obtained difference by a preset proportional coefficient to obtain an initial threshold, traversing all known pollution migration paths to count the ratio of the number of paths with cumulative impedance values ​​exceeding the initial threshold to the total number of known pollution migration paths, and when the ratio exceeds a preset fault tolerance ratio, increasing the preset proportional coefficient by a fixed step size and re-performing the product operation to update the initial threshold, and when the ratio does not exceed the preset fault tolerance ratio, determining the current initial threshold as the preset impedance threshold.

[0014] As a further aspect of the present invention, the specific steps of S5 are as follows: S501: Call the candidate source tracing path, collect the node coordinates within the candidate source tracing path, calculate and extract the path node spacing based on the node coordinates, obtain the preset pollutant attenuation rate, multiply the path node spacing by the preset pollutant attenuation rate to extract attenuation fusion features, and establish a path attenuation index based on the attenuation fusion features. S502: Call the path decay index to obtain the confidence parameter, use the confidence parameter as the numerator and the path decay index as the denominator to perform a ratio operation to extract the decay reciprocal mapping feature, assign a confidence value to the decay reciprocal mapping feature of each independent path unit, perform network weighting on the candidate tracing path based on the confidence value, and establish a path confidence score set. S503: Call the path confidence score set, sort the values ​​in the path confidence score set in descending order and extract the first and second labels, call the corresponding path entity according to the first and second labels, bind the corresponding path entity with business features to extract decision data stream, and generate a source tracing decision scheme based on the decision data stream.

[0015] A knowledge graph-based decision-making system for tracing the source of contaminated sites includes: The water flow component analysis module collects the slope angle, groundwater flow direction angle, and water flow velocity of the contaminated site. Based on the slope angle, it performs spatial direction difference extraction on the groundwater flow direction angle to obtain the flow direction deviation angle. Based on the flow direction deviation angle, it performs direction vector projection analysis on the water flow velocity to generate effective water flow components. The dynamic coupling analysis module calls the effective water flow component to obtain the slope and spatial distance of adjacent land units. Based on the spatial distance, it performs spatial change rate extraction on the slope to obtain the terrain gradient. It then performs dynamic feature coupling processing on the terrain gradient and the effective water flow component to generate a dynamic coupling coefficient. The migration relationship construction module calls the dynamic coupling coefficient, collects the benchmark water level and real-time groundwater level of the contaminated site, performs water level compensation adjustment on the dynamic coupling coefficient, obtains the relationship edge weight, performs feature mapping between adjacent land parcel units based on the relationship edge weight, and constructs pollution migration relationship edges. The graph optimization and comparison module performs network topology splicing processing on the pollution migration relationship edges and adjacent land parcel units to construct a contaminated site knowledge graph. The contaminated site knowledge graph is then input into the Dijkstra algorithm for optimization and comparison to generate candidate source tracing paths. The source tracing decision generation module performs spatial decay feature fusion on the path node spacing of the candidate source tracing paths and the preset pollutant decay rate to obtain the path decay index. Based on the path decay index, it performs confidence assessment and ranking screening on the candidate source tracing paths to generate a source tracing decision scheme.

[0016] Compared with the prior art, the advantages and positive effects of the present invention are as follows: In this invention, the direction vector projection analysis of water flow velocity based on the flow deviation angle generates effective components. Combined with topographic gradient, dynamic coupling processing is performed. Real-time groundwater level compensation relationship edge weights are introduced and topological reconstruction is completed. A knowledge graph of contaminated sites is constructed to drive optimization and comparison. Spatial attenuation features are integrated to obtain the path attenuation index and perform confidence assessment. Through multi-source dynamic feature coupling, the resistance distribution network of pollution migration is reconstructed. The macroscopic map structure and microscopic hydrological dynamic attenuation mechanism are deeply integrated. This breaks away from the operational limitations of single historical comparison and surface interpolation inference, and realizes reliable optimization and ranking of the topological nodes of the full source tracing path, effectively ensuring the rigor and accuracy of the tracing decision results. Attached Figure Description

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

[0018] Figure 1 This is a schematic diagram of the steps of the present invention; Figure 2 This is a detailed schematic diagram of S1 of the present invention; Figure 3 This is a detailed schematic diagram of S2 of the present invention; Figure 4 This is a detailed schematic diagram of S3 of the present invention; Figure 5 This is a detailed schematic diagram of S4 of the present invention; Figure 6 This is a detailed schematic diagram of S5 of the present invention; Figure 7 This is a system module diagram of the present invention. Detailed Implementation

[0019] The technical solution of the present invention will now be described with reference to the accompanying drawings.

[0020] To make the technical problems, technical solutions and advantages of the present invention clearer, a detailed description will be given below in conjunction with the accompanying drawings and specific embodiments.

[0021] Please see Figure 1 This invention provides a knowledge graph-based method for tracing the source of contaminated sites, comprising the following steps: S1: Collect the slope angle, groundwater flow direction angle and water flow velocity of the contaminated site. Based on the slope angle, perform spatial direction difference extraction on the groundwater flow direction angle to obtain the flow direction deviation angle. Based on the flow direction deviation angle, perform direction vector projection analysis on the water flow velocity to generate effective water flow components. S2: Call the effective water flow component, obtain the slope and spatial distance of adjacent plot units, extract the spatial change rate of the slope based on the spatial distance, obtain the terrain gradient, perform dynamic feature coupling processing on the terrain gradient and the effective water flow component, and generate the dynamic coupling coefficient. S3: Call the dynamic coupling coefficient, collect the benchmark water level and real-time groundwater level of the contaminated site, perform water level compensation adjustment on the dynamic coupling coefficient, obtain the relation edge weights, perform feature mapping between adjacent plot units based on the relation edge weights, and construct pollution migration relation edges; S4: Perform network topology splicing on the pollution migration relationship edges and adjacent land parcel units to construct a contaminated site knowledge graph. Input the contaminated site knowledge graph into the Dijkstra algorithm for optimization and comparison to generate candidate source tracing paths. S5: Perform spatial decay feature fusion on the path node spacing of candidate source tracing paths and the preset pollutant decay rate to obtain the path decay index. Based on the path decay index, perform confidence assessment and ranking screening on candidate source tracing paths to generate a source tracing decision scheme.

[0022] The effective water flow component includes the actual migration driving component after the groundwater flow direction is coupled with the slope direction. The dynamic coupling coefficient includes hydraulic gradient, permeability coefficient and hydraulic conductivity. The relation edge weight is specifically the pollution migration association intensity value between plots after water level compensation correction. The contaminated site knowledge graph includes plot nodes, migration relation edges and edge weight distribution. The source tracing decision scheme includes candidate path priority ranking, path confidence score and pollution source pointing result.

[0023] Please see Figure 2 The specific steps of S1 are as follows: S101: Collect the horizontal elevation slope angle, groundwater flow direction angle and instantaneous water flow velocity of the test points of the contaminated site. Extract the projection component of the horizontal elevation slope angle by spatial orthogonal decomposition. Perform Hadamard product operation on the groundwater flow direction angle and the projection component to obtain the site angle distribution matrix. High-precision tilt sensors and Doppler current profilers were deployed at test points in the contaminated site to collect data on horizontal elevation slope angle, groundwater flow direction angle, and instantaneous water flow velocity. During data acquisition, sensor output data was recorded every 10 minutes, and outlier values ​​exceeding the measurement range were removed. Missing data was filled using linear interpolation. The horizontal elevation slope angle refers to the angle between the surface normal and the plumb line at the test point. The groundwater flow direction angle refers to the clockwise angle between the water flow vector projection onto the horizontal plane and true north. The instantaneous water flow velocity was calculated by transmitting sound waves and receiving the echo frequency shift using the Doppler current profiler. During a data acquisition process, the horizontal elevation slope angle at test point 1 was 30.0 degrees, the groundwater flow direction angle was 45.0 degrees, and the instantaneous water flow velocity was 0.5 m / s. At test point 2, the horizontal elevation slope angle was 60.0 degrees, the groundwater flow direction angle was 90.0 degrees, and the instantaneous water flow velocity was 0.8 m / s. After obtaining the horizontal elevation slope angles, they were decomposed spatially orthogonally into a lateral projection component along the east direction and a longitudinal projection component along the north direction. Specifically, the sine of the horizontal elevation slope angle was calculated as the lateral projection component, and the cosine was calculated as the longitudinal projection component. When the horizontal elevation slope angle at test point 1 was 30.0 degrees, the lateral projection component was 0.5, and the longitudinal projection component was 0.866. When the horizontal elevation slope angle at test point 2 was 60.0 degrees, the lateral projection component was 0.866, and the longitudinal projection component was 0.5. After extracting the projected components, the groundwater flow direction angle is multiplied by the corresponding projected component using a Hadamard product operation, which involves multiplying the corresponding elements one by one. For the first test point, multiplying the 45.0-degree groundwater flow direction angle by the lateral projected component (0.5) yields a lateral site angle distribution value of 22.5, and multiplying it by the longitudinal projected component (0.866) yields a longitudinal site angle distribution value of 38.97. For the second test point, multiplying 90.0 degrees by 0.866 yields a lateral site angle distribution value of 77.94, and multiplying it by 0.5 yields a longitudinal site angle distribution value of 45.0. The lateral and longitudinal site angle distribution values ​​of all test points are combined to obtain the site angle distribution matrix.

[0024] S102: Perform discrete partial derivative calculation on adjacent node elements in the site angle distribution matrix to extract feature difference sequence, sort the elements of the feature difference sequence by absolute value comparison, remove abnormal discrete values ​​at both ends with a preset ratio, and perform inverse trigonometric transformation of vector angle on the sorted elements to generate flow direction deviation angle; After obtaining the site angle distribution matrix, discrete partial derivatives are calculated for adjacent node elements within the matrix to extract the feature difference sequence. The logic for calculating discrete partial derivatives is as follows: the difference between adjacent element values ​​is calculated along both the row and column directions of the matrix, and the resulting row difference and column difference are summed to obtain the feature difference. For example, if node element A in the matrix has a value of 22.5, the adjacent right node element B has a value of 23.8, and the adjacent below node element C has a value of 23.8, then the row difference is 23.8 minus 22.5, which equals 1.3, the column difference is 23.8 minus 22.5, which equals 1.3, and the feature difference is 1.3 plus 1.3, which equals 2.6. This calculation is performed on all nodes of the matrix to extract the feature difference sequence. Subsequently, the absolute value of all elements in the feature difference sequence is calculated, and the values ​​are sorted in ascending order. In the sorted sequence, outlier discrete values ​​at both ends are removed. These preset percentages are determined based on the lower limit of the confidence interval for the site's historical hydrological distribution. By reviewing the dispersion statistics for the past six months, the top 5.0% and bottom 5.0% are selected as the preset percentages. Assuming the sorted feature difference sequence contains 100 elements, the first 5 elements and the last 5 elements are directly removed, retaining the middle 90 valid elements. After removing outliers, the retained sorted elements undergo an inverse trigonometric transformation of the vector angle. Specifically, the feature difference is used as the numerator of the spatial location derivative change, and the actual measurement step size of the site's basic grid is used as the denominator. A division operation is performed to obtain the spatial rate of change derivative. Then, the arctangent of this derivative value is used to directly extract the physical deflection angle, thereby accurately generating the flow direction deviation angle. Assuming the extracted feature difference is accurate to 2.598, and the actual measurement step size of the site base grid is set to a fixed 1.5 meters, the derivative of the spatial rate of change obtained by dividing the two is 1.732. Performing arctangent angle conversion on it, the corresponding spatial physical deflection angle is extracted as 60.0 degrees. Thus, the flow direction deviation angle of this calculation node is 60.0 degrees.

[0025] S103: Substitute the flow direction deviation angle into the cosine function to extract the cosine mapping coefficient value, perform scalar dot product calculation on the instantaneous water flow velocity and the cosine mapping coefficient value to extract the projected velocity scalar set, and perform moving average smoothing conversion on the projected velocity scalar set according to the preset sliding window length to generate effective water flow components. After extracting the flow direction deviation angle, the angle is substituted into a cosine function to extract the cosine mapping coefficient value. Specifically, the previously obtained flow direction deviation angle is substituted into the transformation system and the cosine value is calculated. For example, if the flow direction deviation angle calculated in the previous analysis is 60.0 degrees, it is substituted into the cosine function, yielding a cosine mapping coefficient value of 0.5. Subsequently, a scalar dot product is performed on the instantaneous flow velocity and this cosine mapping coefficient value to extract the projected velocity scalar set. The scalar dot product calculation logic is as follows: the instantaneous flow velocity values ​​of each test point are multiplied by the corresponding cosine mapping coefficient value. Assuming the instantaneous flow velocity at a test point is 0.8 m / s and the cosine mapping coefficient value is 0.5, the projected velocity scalar value is 0.8 multiplied by 0.5, which equals 0.4 m / s. The projected velocity scalar values ​​of all test points are summarized to extract the projected velocity scalar set. After extraction, the set is smoothed by a moving average based on a preset sliding window length to generate the effective flow components. The process for determining the preset sliding window length is as follows: The average period of flow velocity fluctuations within the site over the past 30 days is statistically analyzed. The average period duration is divided by 3 to obtain an integer, which is used as the number of window nodes. Here, the window length is set to 5 nodes. The smoothing calculation logic is as follows: Starting from the first position of the set, five consecutive scalar values ​​are sequentially extracted and their arithmetic average is calculated to obtain a smoothed representative value. This smoothed representative value is then calculated sequentially for all positions as the window slides forward. For example, if the first five values ​​in the set are 0.4 m / s, 0.5 m / s, 0.6 m / s, 0.4 m / s, and 0.5 m / s, respectively, their sum is 2.4 m / s. Dividing this by the window length of 5 yields 0.48 m / s, which is taken as the effective flow component at that initial position. As the window slides further, the effective flow components for the entire site are generated sequentially.

[0026] Please see Figure 3 The specific steps of S2 are as follows: S201: Collect the topographic slope and lateral spatial distance of adjacent plots in the contaminated site, perform discrete difference on the topographic slope to extract the node slope difference value sequence, divide the node slope difference value sequence by the lateral spatial distance to extract spatial variation features, and obtain the spatial variation rate matrix. The terrain slope and lateral spatial distance of adjacent plots in the contaminated site were collected using topographic mapping radar equipment. A plane coordinate system was established with the southwest corner of the site as the origin, and measurements were taken sequentially at 5.0-meter intervals along the lateral direction (due east). The terrain slope refers to the inclination of the plot surface, and the lateral spatial distance is the horizontal straight-line length between adjacent test points. After data acquisition, discrete difference analysis was performed on the terrain slope to extract a sequence of node slope difference values. Using the terrain slope of adjacent plots in the contaminated site as input, the slope difference value at each node was calculated by subtracting the slopes of adjacent nodes in the spatial grid. The calculation direction of the slope difference value was sequentially performed along the extension direction of the lateral spatial distance. For example, the terrain slope of plot 1 is 12.0 degrees, plot 2 is 15.0 degrees, and plot 3 is 10.0 degrees, with a lateral spatial distance of 5.0 meters for each plot. The slope difference at node 1 is 15.0 degrees minus 12.0 degrees, equaling 3.0 degrees. The slope difference at node 2 is 10.0 degrees minus 15.0 degrees, equaling -5.0 degrees. The order of elements in the node slope difference value sequence corresponds one-to-one with the order of nodes in the horizontal space of adjacent plots. To clearly illustrate the data distribution rules, Table 1, the test point slope difference value sequence data table, compiles some node calculation data.

[0027] Table 1: Data Table of Slope Difference Values ​​at Test Points

[0028] As shown in Table 1, after extracting the slope difference values ​​for each node, each element in the node slope difference value sequence is divided by a fixed lateral spatial distance to extract spatial variation features and obtain a spatial variation rate matrix. The calculation logic is as follows: the slope difference value of the first node (3.0 degrees) divided by 5.0 meters yields a spatial variation rate of 0.6 degrees per meter; the slope difference value of the second node (-5.0 degrees) divided by 5.0 meters yields a spatial variation rate of -1.0 degrees per meter. This division operation is performed on all difference values ​​in the sequence, and all resulting spatial variation rates are arranged according to their corresponding coordinate positions to generate a spatial variation rate matrix. The advantage of this calculation logic is that by combining discrete difference with distance division, absolute slope is transformed into relative slope gradient, reflecting the degree of dramatic undulation of the terrain in the horizontal direction.

[0029] S202: Perform grid mapping on the spatial rate of change matrix to extract the spatial gradient direction vector, obtain the elevation parameter of the contaminated site, calculate the scalar product of the spatial gradient direction vector and the elevation parameter to extract the initial gradient scalar set, sort the elements in the initial gradient scalar set, and establish the terrain gradient. After extracting the spatial rate of change matrix, a grid mapping is performed on the matrix to extract the spatial gradient direction vector. The grid mapping logic is as follows: each element in the spatial rate of change matrix is ​​mapped to a 2D planar grid model of the site, and the ratio of the rate of change in the longitude and latitude directions of each grid cell is calculated to extract a spatial gradient direction vector with directional attributes. For example, if the spatial rate of change in the longitude direction of a grid cell is 0.6 degrees per meter and the spatial rate of change in the latitude direction is 0.8 degrees per meter, then the length of the spatial gradient direction vector it forms is 1.0 degrees per meter, pointing northeast. Subsequently, the absolute elevation parameter of the contaminated site is obtained through an elevation mapping database, and the above spatial gradient direction vector is scalar-producted with the elevation parameter to extract the initial gradient scalar set. The scalar product calculation logic is as follows: the length of the spatial gradient direction vector is directly multiplied by the elevation parameter value of the corresponding grid node. Assuming the spatial gradient direction vector of a node has a length of 1.0 degrees per meter and the elevation parameter of that node is 120.0 meters, the initial gradient scalar value of that node is calculated by performing a scalar product operation: 1.0 multiplied by 120.0 equals 120.0. This same multiplication operation is performed on all grid nodes to extract an initial gradient scalar set containing the calculation results of all nodes. After extraction, all elements in the initial gradient scalar set are sorted in descending order. The sorting process is performed using the quicksort algorithm, arranging all values ​​from largest to smallest. After sorting, each value and its corresponding original grid node coordinates are strongly correlated and bound, thereby establishing a terrain gradient with a complete spatial hierarchy.

[0030] S203: Based on the effective water flow component and the topographic gradient, perform tensor product operation to extract the dynamic distribution feature matrix, obtain the fluid density constant of the contaminated site, multiply the dynamic distribution feature matrix with the fluid density constant to extract the absolute kinetic energy index set and perform normalization scaling to generate the dynamic coupling coefficient. After obtaining the effective water flow components and topographic gradients, a tensor product operation is performed based on these two parameters. The tensor product operation logic is as follows: the one-dimensional row vector formed by the effective water flow components is multiplied by the one-dimensional column vector formed by the topographic gradients, generating a 2D matrix containing both water flow and topographic dimensions. Assuming that the effective water flow component vector of a certain area contains two elements, 0.48 m / s and 0.52 m / s, and the topographic gradient vector contains two elements, 120.0 and 150.0, after performing the tensor product operation, 0.48 multiplied by 120.0 equals 57.6, 0.48 multiplied by 150.0 equals 72.0, 0.52 multiplied by 120.0 equals 62.4, and 0.52 multiplied by 150.0 equals 78.0. Thus, the four elements of the dynamic distribution feature matrix are extracted as 57.6, 72.0, 62.4, and 78.0. Subsequently, the fluid density constant of the contaminated site was obtained by consulting the site geological survey report. For example, the average fluid density constant of contaminated groundwater was taken as 1.05 kg / L. Each element in the aforementioned dynamic distribution characteristic matrix was multiplied by the fluid density constant to extract the set of absolute kinetic energy indicators. For matrix element 57.6, multiplying by 1.05 yields an absolute kinetic energy indicator of 60.48; for 72.0, multiplying by 1.05 yields 75.6. This process was repeated for all elements to complete the product calculation and combine them into the absolute kinetic energy indicator set. To eliminate differences in data dimensions, normalization scaling was performed on the absolute kinetic energy indicator set. The normalization scaling logic was as follows: find the maximum and minimum values ​​in the set, subtract the minimum value from each indicator, and then divide by the difference between the maximum and minimum values. Assuming the maximum value in the set is 100.0, the minimum value is 50.0, and the current value is 60.48, then the normalized result is 60.48 minus 50.0, which is 10.48. Dividing this by 50.0 gives 0.2096. Outputting all the normalized values ​​generates the final dynamic coupling coefficient.

[0031] Please see Figure 4 The specific steps of S3 are as follows: S301: Collect the benchmark water level and real-time groundwater level of the contaminated site, perform subtraction calculation based on the benchmark water level and real-time groundwater level to extract the drop sequence, perform division calculation based on the drop sequence and benchmark water level to extract the fluctuation ratio set, perform exponential transformation on the fluctuation ratio set, and establish the water level compensation coefficient. Multiple level gauges were deployed at various points within the site to continuously collect baseline water levels and real-time groundwater levels. The baseline water level was defined as the average stable water level height of the site over three consecutive years during the non-high-water season, while the real-time groundwater level was the liquid surface elevation reading directly returned by the sensors within the current detection period. To clearly illustrate the data comparison, Table 2, the water level monitoring characteristic data table, lists the extracted data from some monitoring points.

[0032] Table 2: Water Level Monitoring Characteristic Data Table

[0033] As shown in Table 2, after acquiring the data, a subtraction calculation is performed between the benchmark water level and the real-time groundwater level to extract the head difference sequence. The subtraction calculation logic is as follows: subtract the corresponding benchmark water level value from the real-time groundwater level value to obtain the deviation height of the water level. For example, the real-time groundwater level at monitoring point 1 is 15.5 meters, and the benchmark water level is 15.0 meters; the difference between the two is 0.5 meters. The real-time groundwater level at monitoring point 2 is 14.8 meters, and the benchmark water level is 15.2 meters; the difference between the two is -0.4 meters. The head difference values ​​of the entire site are integrated to extract the head difference sequence. Subsequently, a division calculation is performed between the head difference sequence and the benchmark water level to extract the fluctuation ratio set. The division logic is as follows: divide each calculated head difference value by the benchmark water level value of the corresponding monitoring point. For monitoring point 1, the drop value of 0.5 meters is divided by the baseline water level of 15.0 meters to obtain a fluctuation ratio of 0.033. For monitoring point 2, -0.4 meters is divided by 15.2 meters to obtain a fluctuation ratio of -0.026. These fluctuation ratios are then combined to generate a set. Finally, an exponential transformation is performed on all elements within the fluctuation ratio set to establish water level compensation coefficients. The exponential transformation logic is as follows: using the natural constant as the base and each fluctuation ratio value as the exponent, a power calculation is performed. For a fluctuation ratio of 0.033, the natural constant raised to the power of 0.033 is approximately 1.033; for -0.026, the natural constant raised to the power of -0.026 is approximately 0.974. All calculated power values ​​are output as water level compensation coefficients.

[0034] S302: Multiply the dynamic coupling coefficient and the water level compensation coefficient to extract the connectivity dynamic parameters, perform range standardization mapping on the connectivity dynamic parameters, extract the interval feature values, obtain the connection damping constant, and perform difference operation based on the interval feature values ​​and the connection damping constant to obtain the relation edge weights; After obtaining the dynamic coupling coefficient and water level compensation coefficient output from the aforementioned process, a multiplication operation is performed on these two coefficients to extract the connectivity dynamic parameters. The multiplication logic is as follows: the dynamic coupling coefficient value and the water level compensation coefficient value belonging to the same grid node or plot unit are multiplied one by one. For example, if the dynamic coupling coefficient of a plot unit is calculated to be 0.2096, and its corresponding water level compensation coefficient is 1.033, multiplying the two yields a connectivity dynamic parameter of 0.2165. After traversing all cells in the site to complete the corresponding product calculation, a range standardization mapping is performed on all generated connectivity dynamic parameters to extract interval feature values ​​under a unified dimension. The range standardization mapping logic is as follows: the maximum and minimum values ​​of all connectivity dynamic parameters are retrieved, the minimum value is subtracted from a single connectivity dynamic parameter as the numerator, the difference between the maximum and minimum values ​​is used as the denominator, and the two are divided to obtain the mapped value. Assuming the maximum value among all parameters is 0.8000 and the minimum value is 0.1000, for the connectivity dynamic parameter 0.2165, the numerator is 0.2165 minus 0.1000, equaling 0.1165, and the denominator is 0.8000 minus 0.1000, equaling 0.7000. Dividing these values ​​yields the interval characteristic value of 0.1664. Subsequently, the connection damping constant, which characterizes the resistance of the geological medium to fluid penetration, is obtained by consulting the site soil porosity test documents. Through multi-point sampling and measurement, the average connection damping constant for a typical silty clay layer at this site is set to 0.0500. After obtaining this constant, a difference operation is performed between the calculated interval characteristic value and the connection damping constant to obtain the final relation edge weights. The difference operation logic is to directly subtract the connection damping constant from the interval characteristic value. For the interval characteristic value of 0.1664, subtracting the connection damping constant of 0.0500, the weight of the relation edge connecting the cell node to the outside is 0.1164.

[0035] S303: Call the relation edge weights, extract the spatial topological coordinate matrix of adjacent plot units, concatenate the relation edge weights with the topological coordinate matrix to extract the associated feature vectors, sort the associated feature vectors in descending order and remove the end nodes, assign the remaining connections to the corresponding parameters, and generate pollution migration relation edges. The system calls upon the relational edge weight dataset and simultaneously extracts the spatial topological coordinate matrix of adjacent land parcel units from the spatial geographic information database. This spatial topological coordinate matrix completely records the 3D coordinates of the centers of all land parcels divided into independent units within the contaminated site, namely longitude, latitude, and elevation feature data. For example, the topological coordinates of a certain central land parcel are longitude 116.0 degrees, latitude 39.0 degrees, and elevation 15.0 meters. After obtaining the coordinate matrix, the previously calculated relational edge weights are concatenated with this topological coordinate matrix to extract the associated feature vector. The concatenation logic is as follows: the relational edge weight value corresponding to each group of adjacent land parcels is used as the fourth dimension data and appended to the start and end coordinate information of that group of adjacent land parcels to form a complete directed feature record. Assuming that the relational edge weight from the starting land parcel to the target land parcel is 0.1164, it is merged with the 3D coordinates of the two land parcels to form a single associated feature vector. After completing the concatenation and combination of the entire site, all generated associated feature vectors are sorted in descending order according to the magnitude of the relational edge weight values. In a long sequence arranged in descending order, the connection data corresponding to the terminal nodes are removed. The removal ratio of terminal nodes is set to the last 10.0% of the sorted sequence, based on the fact that the probability of ineffective connectivity through weak pores measured in historical penetration experiments is approximately 9.5%, hence the integer 10.0% is used as the truncation threshold. Assuming a total of 1000 associated feature vectors are generated, after sorting them in descending order of weight, the last 100 extremely weak associated vectors are directly deleted. After the removal process, the remaining 900 connection data are assigned corresponding coordinates and weight parameters, formally generating contamination migration relationship edges in the data space.

[0036] Please see Figure 5 The specific steps of S4 are as follows: S401: Call the pollution migration relationship edges, collect spatial topological feature parameters, perform connection reconstruction based on the spatial topological feature parameters to extract the topological connection matrix, define node entities according to the topological connection matrix, aggregate the node entities with the pollution migration relationship edges, and establish a knowledge graph of the contaminated site. After extracting the pollution migration relationship edges, the dataset of these edges is retrieved, and spatial topological feature parameters are collected through the site environmental monitoring platform. These parameters encompass the spatial distribution information of various geographical obstacles and artificial isolation structures, such as the location of anti-seepage walls, bedrock uplift areas, and the direction of underground centralized pipelines. Based on these spatial topological feature parameters, a connection reconstruction operation is performed on the existing connected network to extract a new topological connection matrix. The connection reconstruction logic involves checking each pollution migration relationship edge for spatial intersection with the physical blocking elements identified by the spatial topological feature parameters. If a spatial intersection occurs, the connectivity of that migration relationship edge is reset to disconnected. For example, if the starting and ending points of a pollution migration relationship edge cross an underground anti-seepage wall 20.0 meters deep, the adjacency value of this edge in the topological connection matrix is ​​forcibly changed from 1 to 0, indicating physical disconnection. Edges not blocked retain an adjacency value of 1. After full-domain data comparison and state reset, a topological connection matrix reflecting the true geological connectivity state is output. Based on the row and column distribution characteristics of the topological connection matrix, basic node entities in the network are defined. Each node entity not only includes basic three-dimensional coordinates but also incorporates static characteristic attribute data such as soil texture parameters and historical sampling concentrations at that location. After defining the multidimensional attributes of the node entities, a data aggregation operation is performed between the scattered node entities and the previously reconstructed and verified pollution migration relationship edges. The aggregation process involves using the node entities as vertices of the graph in the graph database, and the pollution migration relationship edges and their associated weights as directed edges, utilizing the mapping relationship between vertices and edges to load the entire network data. This ultimately constructs a well-structured knowledge graph of contaminated sites with constraints based on physical entity meaning.

[0037] S402: Call the contaminated site knowledge graph, read historical test data to extract contaminated nodes and target nodes, retrieve the connection paths between nodes according to the topology connection matrix, accumulate the contaminated migration relationship edge association parameters step by step to extract the cumulative impedance value, and establish a path impedance weight set. After establishing a knowledge graph of the contaminated site, this graph is invoked and historical test data is read to extract confirmed contamination nodes with high concentrations exceeding standards, as well as downstream target nodes awaiting investigation. For example, historical underground test data shows that the chemical oxygen demand (COD) of sampling well No. 3 on the north side of the site is severely exceeded, so it is directly marked as a contamination node, and monitoring wells No. 8 and No. 9 on the south boundary of the site are marked as the final target nodes. Based on the topological connectivity matrix embedded in the graph, all potential spatial connectivity paths from the aforementioned contamination nodes to the target nodes are retrieved from the graph database. A connectivity path refers to a continuous flow channel in the graph composed of a series of directed edges connected end to end. For each retrieved complete path, the correlation parameters of the contamination migration relationship edges constituting the path are accumulated level by level to extract the global cumulative impedance value of the path. The correlation parameter is the reciprocal of the previously set relation edge weights; the larger the weight of a single segment, the easier the material migration, and the smaller the corresponding impedance. The accumulation logic is to sum the impedance values ​​of each directed edge segment along the spatial path starting from the contamination node. Suppose a spatially connected path consists of three directed edges. The first edge has a weight of 0.2 and an impedance value of the reciprocal of its weight, 5.0. The second edge has a weight of 0.4 and an impedance value of 2.5. The third edge has a weight of 0.5 and an impedance value of 2.0. The cumulative impedance of this path is 5.0 + 2.5 + 2.0 = 9.5. By iterating through all retrieved connected paths and repeating the above step-by-step accumulation operation, the total impedance of each path is calculated. To clearly illustrate the distribution of path weights, Table 3, the path impedance weight set table, lists the calculation results for some channels.

[0038] Table 3: Set of Path Impedance Weights

[0039] As shown in Table 3, the unique numbers corresponding to all paths and the final calculated cumulative impedance values ​​are packaged and compiled to establish a set of path impedance weights with global comparison and optimization significance.

[0040] S403: Call the path impedance weight set, sort the values ​​in ascending order and extract the first value, define the path corresponding to the first value as the benchmark connected path, compare the remaining values ​​with the preset impedance threshold to remove items exceeding the threshold, obtain the set of qualified connected paths, and concatenate the paths in the set of qualified connected paths with the benchmark connected paths to extract the optimized related subset and generate candidate traceability paths. The path impedance weight set in Table 3 is called, and all cumulative impedance values ​​in the set are sorted in ascending order. After sorting, the first value at the top is extracted. Since the lowest impedance represents the smoothest material migration, the physical path corresponding to the first value is directly defined as the baseline connected path. Assuming the first value in the above sorting is 8.2, then its corresponding third path is the baseline connected path. Subsequently, the remaining cumulative impedance values ​​in the set are compared with a preset impedance threshold to eliminate invalid path items that exceed the threshold limit. The process of setting the preset impedance threshold is as follows: First, the cumulative impedance values ​​corresponding to all known pollution migration paths in the historical test data are extracted, and the maximum and minimum values ​​are calculated. Assuming the historical maximum impedance is 20.0 and the minimum impedance is 5.0, the difference is 15.0. The difference 15.0 is multiplied by the preset scaling factor 0.6 to obtain an initial threshold of 9.0. Then, all known pollution migration paths are traversed, and the ratio of the number of paths with impedance values ​​exceeding the initial threshold of 9.0 to the total number of paths is counted. Assume there are a total of 50 known paths, of which 10 exceed 9.0, a ratio of 20.0%. When this ratio exceeds the preset tolerance ratio of 15.0%, the preset ratio coefficient of 0.6 is increased to 0.65 in fixed steps of 0.05, and then multiplied by the difference of 15.0 to update the initial threshold to 9.75. The excess ratio is then counted again. If the number of excess paths decreases to 6, with a ratio of 12.0%, no longer exceeding the preset tolerance ratio of 15.0%, then the iterative update stops, and the current initial threshold of 9.75 is determined as the final preset impedance threshold. The remaining paths are filtered using the preset impedance threshold of 9.75, for example, removing the second path with an impedance of 12.4 and retaining the first path with an impedance of 9.5. All retained paths less than or equal to 9.75 are extracted to obtain the set of compliant connected paths. Finally, all paths in the qualified connectivity path set are concatenated and recombined with the initially determined baseline connectivity path to extract the optimization correlation subset containing the primary and secondary preferred channels, and finally output as candidate tracing paths for substantive tracing decisions.

[0041] Please see Figure 6 The specific steps of S5 are as follows: S501: Call the candidate source tracing path, collect the node coordinates within the candidate source tracing path, calculate and extract the path node spacing based on the node coordinates, obtain the preset pollutant attenuation rate, multiply the path node spacing with the preset pollutant attenuation rate to extract attenuation fusion features, and establish a path attenuation index based on the attenuation fusion features. After extracting candidate source tracing paths, the dataset of these paths is accessed, and the 3D coordinates of all nodes within each path are collected from all angles. Based on the 3D spatial coordinates of adjacent nodes, the straight-line distance between nodes in the physical path, i.e., the path node spacing, is calculated using the Euclidean distance principle. The calculation logic is as follows: the coordinate differences between two adjacent nodes in the three dimensions of longitude, latitude, and elevation are taken, their squares are calculated, and then summed. Finally, the square root is taken. For example, assuming the coordinates of the first node are 100.0 meters longitude, 100.0 meters latitude, and 10.0 meters elevation, and the coordinates of the second node are 103.0 meters longitude, 104.0 meters latitude, and 10.0 meters elevation, the longitude and latitude differences are 3.0 meters and 4.0 meters respectively. The sum of their squares is 9.0 plus 16.0, which equals 25.0. Taking the square root yields a path node spacing of 5.0 meters. Subsequently, the preset pollutant attenuation rate is obtained by reviewing the site's previously measured underground pollutant degradation cycle report. The preset pollutant attenuation rate characterizes the ratio of natural degradation to physical adsorption loss of a specific pollutant per unit distance migrated in the underground stratum medium. For the typical heavy metal chromium pollutant in this contaminated site, the preset pollutant attenuation rate was determined to be 0.2% per meter based on previous experiments. After obtaining this preset pollutant attenuation rate, it is multiplied by the previously calculated path node spacing to extract the attenuation fusion feature. The multiplication logic is as follows: multiply the single path node spacing of 5.0 meters by the preset pollutant attenuation rate of 0.2% per meter to calculate the attenuation fusion feature value between nodes in that segment as 1.0%. After performing the same product calculation on each adjacent node in the candidate source tracing path, the attenuation fusion feature values ​​of each segment on the entire path are multiplied together or summed to serve as a macroscopic damping metric. Based on this macroscopic metric, a path attenuation index with a specific dimension is established for each candidate path.

[0042] S502: Call the path decay index to obtain the confidence parameter. Use the confidence parameter as the numerator and the path decay index as the denominator to perform a ratio operation to extract the decay inverse mapping feature. Assign a confidence value to each independent path unit based on the decay inverse mapping feature. Perform network weighting on the candidate source tracing path based on the confidence value to establish a path confidence score set. After generating the path decay indices for all candidates, the path decay index values ​​for each business entry are retrieved. Simultaneously, a global confidence parameter for the specific site evaluation task is obtained from the backend rule base. The initial setting of the confidence parameter is based on the median signal-to-noise ratio (SNR) of independent verification point samples collected during the preliminary soil exploration of the site. After statistically analyzing historical characteristic data from multiple site explorations, a SNR of 10.0 is selected as the global confidence parameter for this core calculation. After obtaining the confidence parameter and path decay indices, the fixed confidence parameter is used as the numerator, and the path decay index corresponding to each candidate path is used as the denominator. A ratio calculation is then performed to extract the inverse decay mapping feature. The ratio calculation logic is to directly divide the numerator by the denominator. For example, assuming that for a candidate tracing path, the aforementioned derivation steps calculate its path decay index to be 2.5, then using the global confidence parameter of 10.0 as the numerator and the path decay index of 2.5 as the denominator, the division yields the inverse decay mapping feature of this path as 4.0. Similarly, if the decay index of another spatial path is 5.0, then 10.0 divided by 5.0 yields a decay inverse mapping feature of 2.0. After the division operation, for each independently calculated path unit, the extracted decay inverse mapping feature value is directly used as a quantifiable scoring evaluation index, mapped and assigned to the corresponding path unit to form the final confidence value. For the first path mentioned above, a confidence value of 4.0 is assigned. After obtaining the confidence values ​​of all independent path units, a network weighting operation is performed on all candidate source tracing paths based on these confidence values. The weighting operation logic is as follows: a normalization function is used to scale all confidence values, limiting their sum to 100.0%, and then the scaled percentage relative weights are bound to the corresponding path edges in the knowledge graph. Assuming that only the above two candidate paths exist in the computational library, with confidence values ​​of 4.0 and 2.0 respectively, and a sum of 6.0, the weight of the first path after scaling is 66.7%, and that of the second path is 33.3%. After completing all weighting actions, the underlying data structure containing the final weight parameters and path mapping relationships will be output in a unified manner to establish a set of path confidence scores for the entire network.

[0043] S503: Call the path confidence score set, sort the values ​​in the path confidence score set in descending order and extract the first and second labels, call the corresponding path entity according to the first and second labels, bind the corresponding path entity with business features to extract the decision data stream, and generate a source tracing decision scheme based on the decision data stream; After establishing the path confidence score set for the entire network, the score weight data of all detailed records within this dataset is retrieved. All network weight score values ​​recorded in the path confidence score set are sorted in descending order. The sorting logic is as follows: the backend uses a fast comparison sorting algorithm to strictly adjust the data positions of the underlying key-value pair structure containing path labels and score values ​​according to the score values ​​from highest to lowest. After the descending sort is completed, the optimal label feature information at the top of the sequence is directly extracted. Assuming the set contains 5 candidate path scores, with the highest score being 66.7%, and its corresponding path feature label being path number 1, this label is precisely extracted. Based on the extracted first label, the complete path entity data pointed to by this label is retrieved from the underlying graph database. The retrieved corresponding path entity contains detailed information including the three-dimensional coordinates, flow direction, medium damping, and the associated weights of all directed edges along the path from the suspected pollution outbreak point to the current monitoring target node. Subsequently, the extracted corresponding path entity is bound to the macro-business features of the site pollution remediation management system. The business characteristics include the land planning nature of different areas within the contaminated site, the upper limit of the remediation budget, and the identification of surrounding sensitive ecological protection areas. The binding logic involves spatially overlaying and comparing each segment of the spatial grid traversed by the physical network path with the 3D model of the site planning, automatically marking the special areas encountered along the path and the required remediation difficulty index. After completing the deep binding of multi-dimensional attribute information, all spatial path feature data and business constraints are encapsulated and packaged to form a highly standardized terminal decision data stream. Finally, a complete source tracing decision plan is generated based on the comprehensive information within this decision data stream. This source tracing decision plan can intuitively and scientifically indicate the exact coordinates of the pollution source and output a profile map of the main channels for pollution migration and an early warning assessment level for the pollution at each sensitive node along the route.

[0044] Please see Figure 7 A knowledge graph-based decision-making system for tracing the source of contaminated sites includes: The water flow component analysis module collects the slope angle, groundwater flow direction angle, and water flow velocity of the contaminated site. Based on the slope angle, it performs spatial direction difference extraction on the groundwater flow direction angle to obtain the flow direction deviation angle. Based on the flow direction deviation angle, it performs direction vector projection analysis on the water flow velocity to generate effective water flow components. The dynamic coupling analysis module calls the effective water flow component to obtain the slope and spatial distance of adjacent plot units. Based on the spatial distance, it performs spatial change rate extraction on the slope to obtain the terrain gradient. It then performs dynamic feature coupling processing on the terrain gradient and the effective water flow component to generate the dynamic coupling coefficient. The migration relationship construction module calls the dynamic coupling coefficient, collects the benchmark water level and real-time groundwater level of the contaminated site, performs water level compensation adjustment on the dynamic coupling coefficient, obtains the relationship edge weight, performs feature mapping between adjacent plot units based on the relationship edge weight, and constructs pollution migration relationship edges. The graph optimization and comparison module performs network topology splicing processing on pollution migration relationship edges and adjacent land parcel units to construct a contaminated site knowledge graph. The contaminated site knowledge graph is input into the Dijkstra algorithm for optimization and comparison to generate candidate source tracing paths. The source tracing decision generation module performs spatial decay feature fusion on the path node spacing of candidate source tracing paths and the preset pollutant decay rate to obtain the path decay index. Based on the path decay index, it performs confidence assessment and ranking screening on candidate source tracing paths to generate source tracing decision schemes.

[0045] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of protection of the described technical solutions.

Claims

1. A pollution site tracing method based on a knowledge graph, characterized in that, Includes the following steps: S1: Collect the slope angle, groundwater flow direction angle and water flow velocity of the contaminated site. Based on the slope angle, perform spatial direction difference extraction on the groundwater flow direction angle to obtain the flow direction deviation angle. Based on the flow direction deviation angle, perform direction vector projection analysis on the water flow velocity to generate effective water flow components. S2: Call the effective water flow component to obtain the slope and spatial distance of adjacent land units, extract the spatial change rate of the slope based on the spatial distance, obtain the terrain gradient, perform dynamic feature coupling processing on the terrain gradient and the effective water flow component, and generate a dynamic coupling coefficient. S3: Call the dynamic coupling coefficient, collect the benchmark water level and real-time groundwater level of the contaminated site, perform water level compensation adjustment on the dynamic coupling coefficient, obtain the relation edge weights, perform feature mapping between adjacent land parcel units based on the relation edge weights, and construct pollution migration relation edges; S4: Perform network topology splicing processing on the pollution migration relationship edge and adjacent land parcel units to construct a contaminated site knowledge graph. Input the contaminated site knowledge graph into the Dijkstra algorithm for optimization and comparison to generate candidate source tracing paths. S5: Perform spatial decay feature fusion on the path node spacing of the candidate source tracing path and the preset pollutant decay rate to obtain the path decay index. Based on the path decay index, perform confidence assessment and ranking screening on the candidate source tracing path to generate a source tracing decision scheme.

2. The knowledge graph-based method for pollution site provenance according to claim 1, wherein, The effective water flow component includes the actual migration driving component after coupling the groundwater flow direction with the slope direction. The dynamic coupling coefficient includes hydraulic gradient, permeability coefficient and hydraulic conductivity. The relation edge weight is specifically the pollution migration association intensity value between plots after water level compensation correction. The contaminated site knowledge graph includes plot nodes, migration relation edges and edge weight distribution. The source tracing decision scheme includes candidate path priority ranking, path confidence score and pollution source pointing result.

3. The knowledge graph-based method for tracing the source of contaminated sites according to claim 1, characterized in that, The specific steps of S1 are as follows: S101: Collect the horizontal elevation slope angle, groundwater flow direction angle and instantaneous water flow velocity of the test points of the contaminated site. Extract the projection component of the horizontal elevation slope angle by spatial orthogonal decomposition. Perform Hadamard product operation on the groundwater flow direction angle and the projection component to obtain the site angle distribution matrix. S102: Perform discrete partial derivative calculation on adjacent node elements in the site angle distribution matrix to extract feature difference sequence, sort the elements of the feature difference sequence by absolute value comparison, remove abnormal discrete values ​​at both ends with a preset ratio, and perform inverse trigonometric transformation of vector angle on the sorted elements to generate flow direction deviation angle; S103: Substitute the flow direction deviation angle into the cosine function to extract the cosine mapping coefficient value, perform scalar dot product calculation on the instantaneous water flow velocity and the cosine mapping coefficient value to extract the projected velocity scalar set, and perform sliding average smoothing conversion on the projected velocity scalar set according to the preset sliding window length to generate effective water flow components.

4. The knowledge graph-based method for tracing the source of contaminated sites according to claim 3, characterized in that, The specific steps of S2 are as follows: S201: Collect the topographic slope and lateral spatial distance of adjacent plots in the contaminated site, perform discrete difference on the topographic slope to extract the node slope difference value sequence, divide the node slope difference value sequence by the lateral spatial distance to extract spatial variation features, and obtain the spatial variation rate matrix. S202: Perform grid mapping on the spatial rate of change matrix to extract the spatial gradient direction vector, obtain the elevation parameter of the contaminated site, perform scalar product calculation on the spatial gradient direction vector and the elevation parameter to extract the initial gradient scalar set, sort the elements in the initial gradient scalar set, and establish the terrain gradient. S203: Based on the effective water flow component and the terrain gradient, perform tensor product operation to extract the dynamic distribution feature matrix, obtain the fluid density constant of the contaminated site, multiply the dynamic distribution feature matrix with the fluid density constant to extract the absolute kinetic energy index set and perform normalization scaling to generate the dynamic coupling coefficient.

5. The knowledge graph-based method for tracing the source of contaminated sites according to claim 4, characterized in that, The discrete difference extraction of node slope difference value sequence for terrain slope refers to taking the terrain slope of adjacent plots in the contaminated site as input, and calculating the slope difference value at each node in the spatial grid by subtracting the slope of adjacent nodes. The calculation direction of the slope difference value is executed sequentially along the extension direction of the horizontal spatial distance, and the arrangement order of each element in the node slope difference value sequence corresponds one-to-one with the arrangement order of nodes of adjacent plots in the horizontal space.

6. The knowledge graph-based method for tracing the source of contaminated sites according to claim 4, characterized in that, The specific steps for S3 are as follows: S301: Collect the benchmark water level and real-time groundwater level of the contaminated site, perform subtraction calculation based on the benchmark water level and real-time groundwater level to extract the drop sequence, perform division calculation based on the drop sequence and benchmark water level to extract the fluctuation ratio set, perform exponential transformation on the fluctuation ratio set, and establish the water level compensation coefficient. S302: Multiply the dynamic coupling coefficient and the water level compensation coefficient to extract the connectivity dynamic parameter, perform range standardization mapping on the connectivity dynamic parameter, extract the interval feature value, obtain the connection damping constant, and perform difference operation based on the interval feature value and the connection damping constant to obtain the relation edge weight; S303: Call the relation edge weights, extract the spatial topological coordinate matrix of adjacent plot units, concatenate the relation edge weights with the topological coordinate matrix to extract the associated feature vectors, sort the associated feature vectors in descending order and remove the end nodes, assign the remaining connections with corresponding parameters, and generate pollution migration relation edges.

7. The knowledge graph-based method for tracing the source of contaminated sites according to claim 6, characterized in that, The specific steps of S4 are as follows: S401: Call the pollution migration relationship edge, collect spatial topological feature parameters, perform connection reconstruction based on the spatial topological feature parameters to extract the topological connection matrix, define node entities according to the topological connection matrix, aggregate the node entities with the pollution migration relationship edge, and establish a knowledge graph of the contaminated site. S402: Call the contaminated site knowledge graph, read historical test data to extract contaminated nodes and target nodes, retrieve the connection paths between nodes according to the topology connection matrix, accumulate the contaminated migration relationship edge association parameters step by step to extract the cumulative impedance value, and establish a path impedance weight set. S403: Call the path impedance weight set, sort the values ​​in ascending order and extract the first value, define the path corresponding to the first value as the benchmark connected path, compare the remaining values ​​with the preset impedance threshold to remove items exceeding the threshold, obtain the set of qualified connected paths, and concatenate the paths in the set of qualified connected paths with the benchmark connected paths to extract the optimized related subset and generate candidate traceability paths.

8. The knowledge graph-based method for tracing the source of contaminated sites according to claim 7, characterized in that, The preset impedance threshold is obtained by extracting the cumulative impedance values ​​corresponding to all known pollution migration paths in historical test data, performing a difference operation between the maximum and minimum values ​​of all cumulative impedance values, multiplying the resulting difference by a preset proportional coefficient to obtain the initial threshold, traversing all known pollution migration paths to count the ratio of the number of paths with cumulative impedance values ​​exceeding the initial threshold to the total number of known pollution migration paths, and when the ratio exceeds a preset fault tolerance ratio, increasing the preset proportional coefficient by a fixed step size and re-performing the product operation to update the initial threshold, and when the ratio does not exceed the preset fault tolerance ratio, the current initial threshold is determined as the preset impedance threshold.

9. The knowledge graph-based method for tracing the source of contaminated sites according to claim 7, characterized in that, The specific steps of S5 are as follows: S501: Call the candidate source tracing path, collect the node coordinates within the candidate source tracing path, calculate and extract the path node spacing based on the node coordinates, obtain the preset pollutant attenuation rate, multiply the path node spacing by the preset pollutant attenuation rate to extract attenuation fusion features, and establish a path attenuation index based on the attenuation fusion features. S502: Call the path decay index to obtain the confidence parameter, use the confidence parameter as the numerator and the path decay index as the denominator to perform a ratio operation to extract the decay reciprocal mapping feature, assign a confidence value to the decay reciprocal mapping feature of each independent path unit, perform network weighting on the candidate tracing path based on the confidence value, and establish a path confidence score set. S503: Call the path confidence score set, sort the values ​​in the path confidence score set in descending order and extract the first and second labels, call the corresponding path entity according to the first and second labels, bind the corresponding path entity with business features to extract decision data stream, and generate a source tracing decision scheme based on the decision data stream.

10. A knowledge graph-based decision-making system for tracing the source of contaminated sites, characterized in that, The system is used to implement the knowledge graph-based source tracing method for contaminated sites as described in any one of claims 1-9, and the system comprises: The water flow component analysis module collects the slope angle, groundwater flow direction angle, and water flow velocity of the contaminated site. Based on the slope angle, it performs spatial direction difference extraction on the groundwater flow direction angle to obtain the flow direction deviation angle. Based on the flow direction deviation angle, it performs direction vector projection analysis on the water flow velocity to generate effective water flow components. The dynamic coupling analysis module calls the effective water flow component to obtain the slope and spatial distance of adjacent land units. Based on the spatial distance, it performs spatial change rate extraction on the slope to obtain the terrain gradient. It then performs dynamic feature coupling processing on the terrain gradient and the effective water flow component to generate a dynamic coupling coefficient. The migration relationship construction module calls the dynamic coupling coefficient, collects the benchmark water level and real-time groundwater level of the contaminated site, performs water level compensation adjustment on the dynamic coupling coefficient, obtains the relationship edge weight, performs feature mapping between adjacent land parcel units based on the relationship edge weight, and constructs pollution migration relationship edges. The graph optimization and comparison module performs network topology splicing processing on the pollution migration relationship edges and adjacent land parcel units to construct a contaminated site knowledge graph. The contaminated site knowledge graph is then input into the Dijkstra algorithm for optimization and comparison to generate candidate source tracing paths. The source tracing decision generation module performs spatial decay feature fusion on the path node spacing of the candidate source tracing paths and the preset pollutant decay rate to obtain the path decay index. Based on the path decay index, it performs confidence assessment and ranking screening on the candidate source tracing paths to generate a source tracing decision scheme.

Images