A method and system for tracing groundwater organic pollution by fusing multi-source heterogeneous data

Abstract

The application relates to the technical field of underground water pollution monitoring, and discloses a method and system for tracing the source of underground water organic pollution by fusing multi-source heterogeneous data. The method comprises the following steps: obtaining attribute fluctuation characteristics by processing geological difference data through a random forest; judging whether space-time dislocation exists in monitoring data fragments according to the attribute fluctuation characteristics; if yes, a first monitoring sequence is obtained by reconstruction processing; if not, concentration values are extracted as a second monitoring sequence; a concentration distribution matrix is obtained by carrying out grid division and spatial interpolation on the monitoring sequence; diffusion evolution processing is carried out on the concentration distribution matrix in combination with the attribute fluctuation characteristics to obtain a spreading track; and a pollution source is obtained by performing reverse tracing on the diffusion path extracted from the spreading track. The method can solve the problem that the fusion of multi-source heterogeneous geological data and monitoring data is insufficient, resulting in low tracing accuracy of the source of underground water organic pollution.

Description

Technical Field

[0001] This invention relates to the field of groundwater pollution monitoring technology, and in particular to a method and system for tracing the source of organic pollution in groundwater by integrating multi-source heterogeneous data. Background Technology

[0002] Tracing the source of organic pollution in groundwater is a core technological step in ensuring ecological and environmental safety and public drinking water health. It relies on data acquisition and monitoring and control technologies to achieve continuous perception and monitoring of the groundwater environment.

[0003] Existing technologies collect time-series data on water level and quality in target areas by deploying sensor monitoring networks, and then simulate the underground migration process of pollutants using hydrogeological parameter models. However, in areas with complex three-dimensional geological structures, the permeability and porosity of different strata exhibit dramatic spatial differences, leading to fragmentation and spatiotemporal misalignment of sensor monitoring data. This results in a lack of effective data linkage between forward migration simulation and reverse source tracing, and a significant deviation between the pollution source location results and the actual source.

[0004] In summary, existing technologies suffer from the problem of low accuracy in tracing the source of organic pollution in groundwater due to the difficulty in effectively integrating three-dimensional geological spatial difference data, hydrogeological attribute fluctuation characteristics, and fragmented monitoring data. Summary of the Invention

[0005] This invention provides a method and system for tracing the source of organic pollution in groundwater by integrating multi-source heterogeneous data, in order to solve the problem of low accuracy in tracing the source of organic pollution in groundwater due to insufficient integration of multi-source heterogeneous data in the prior art.

[0006] Firstly, to address the aforementioned technical problems, this invention provides a method for tracing the source of organic pollution in groundwater by integrating multi-source heterogeneous data, comprising: Geological difference data of the target area is obtained, and random forest processing is applied to the geological difference data to obtain attribute fluctuation characteristics; Acquire monitoring data fragments of the target area, determine whether the monitoring data fragments have spatiotemporal misalignment based on the attribute fluctuation characteristics, if so, reconstruct the monitoring data fragments to obtain a first monitoring sequence, if not, extract the concentration values ​​of the monitoring data fragments as a second monitoring sequence; The first monitoring sequence or the second monitoring sequence is divided into grids and spatially interpolated to obtain the concentration distribution matrix of pollutants. The concentration distribution matrix is ​​then subjected to diffusion evolution processing based on the attribute fluctuation characteristics to obtain the pollutant spread trajectory. Based on the pollutant spread trajectory, the pollutant diffusion path is extracted, and the pollutant diffusion path is traced backward to obtain the pollution source.

[0007] Secondly, the present invention provides a groundwater organic pollution tracing system that integrates multi-source heterogeneous data, comprising: a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements the steps of the groundwater organic pollution tracing method that integrates multi-source heterogeneous data as described above.

[0008] Thirdly, the present invention provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the groundwater organic pollution tracing method as described above, which integrates multi-source heterogeneous data.

[0009] Compared with the prior art, the present invention has the following beneficial effects: (1) This invention extracts attribute fluctuation features from geological difference data through random forest, and reconstructs spatiotemporally misaligned monitoring data based on the attribute fluctuation features, organically linking geological feature analysis and monitoring data repair, solving the problem that geological spatial difference data and fragmented monitoring data are difficult to effectively integrate in the prior art, and realizing the continuity and integrity of monitoring sequence driven by multi-source heterogeneous data.

[0010] (2) This invention corrects the concentration distribution matrix by topographic resistance and diffusion evolution by attribute fluctuation characteristics, introduces geological heterogeneity constraints into the pollutant migration simulation process, so that the generated spread trajectory can truly reflect the differentiated blocking effect of strata with different permeability on pollutant diffusion, and improves the spatial accuracy and physical consistency of forward spread simulation.

[0011] (3) This invention identifies bifurcation nodes by calculating the connectivity and deflection angle of the pollutant diffusion path, and implements differentiated reverse tracing strategies for the two scenarios of having bifurcation and not having bifurcation, thereby avoiding the problem of the tracing path extending along the wrong branch due to underground bypass flow and improving the accuracy of pollution source location. Attached Figure Description

[0012] Figure 1 This is a flowchart illustrating a method for tracing the source of organic pollution in groundwater by integrating multi-source heterogeneous data, as provided in an embodiment of the present invention. Detailed Implementation

[0013] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0014] Reference Figure 1The first embodiment of the present invention provides a method for tracing the source of organic pollution in groundwater by fusing multi-source heterogeneous data, including the following steps: S11, Obtain geological difference data of the target area, and perform random forest processing on the geological difference data to obtain attribute fluctuation characteristics; S12, acquire monitoring data fragments of the target area, determine whether the monitoring data fragments have spatiotemporal misalignment based on the attribute fluctuation characteristics, if so, reconstruct the monitoring data fragments to obtain a first monitoring sequence, if not, extract the concentration values ​​of the monitoring data fragments as a second monitoring sequence. S13, perform grid division and spatial interpolation on the first monitoring sequence or the second monitoring sequence to obtain the concentration distribution matrix of pollutants, and perform diffusion evolution processing on the concentration distribution matrix according to the attribute fluctuation characteristics to obtain the pollutant spread trajectory; S14. Extract the pollutant diffusion path based on the pollutant spread trajectory, and perform reverse tracing on the pollutant diffusion path to obtain the pollution source.

[0015] In step S11, geological difference data of the target area is acquired, and random forest processing is applied to the geological difference data to obtain attribute fluctuation characteristics, including: The geological difference data are subjected to feature extraction and normalization to obtain the input feature set; Random forest regression prediction is performed on the input feature set to obtain the water level change sequence; The amplitude and frequency of the water level change sequence are extracted to obtain the attribute fluctuation characteristics.

[0016] Specifically, through a comprehensive exploration approach including borehole drilling, ground-penetrating radar (GPR), and geophysical logging, raw geological parameters such as stratum thickness, porosity, and permeability are collected at each exploration point within the target area. The stratum thickness data is used to calculate the elevation difference (vertical distance between adjacent strata interfaces) of each grid node and to determine the aquifer layer structure, while simultaneously obtaining the thickness data of the target aquifer. In this embodiment, the average thickness of the unconfined aquifer in the target area is approximately 1 m, ranging from 0.8 m to 1.3 m. This thickness is obtained by cross-checking the borehole columnar section and geophysical logging interpretation results, and after processing, it is incorporated into the subsequent PCA input feature set and flow tube equivalent diameter setting in the form of elevation difference and aquifer layer thickness.

[0017] It should be noted that the discrete borehole data is spatially gridded using the Kriging interpolation algorithm, with a grid resolution of 50m × 50m. The porosity gradient and permeability difference are calculated at each grid node. The porosity gradient is the rate of change of porosity between that node and its neighboring nodes, and the permeability difference is the absolute difference in permeability between that node and its neighboring nodes. These porosity gradient and permeability difference values ​​are then combined to output the geological difference data for the corresponding grid node.

[0018] For example, two adjacent grid nodes A and B, 100m apart horizontally in the target aquifer, are selected for illustration. The permeability at node A is 12.5 m / d, and the porosity is 0.28; the permeability at node B is 11.7 m / d, and the porosity is 0.25. Calculations show that the permeability variation rate in this direction is 0.8 m / d per 100m, and the porosity variation rate is 0.03 per 100m, constituting the geological difference data at this location.

[0019] Specifically, the permeability change rate (X / Y / Z components), porosity change rate (X / Y / Z components), elevation difference, and historical meteorological and hydrological characteristics (including daily rainfall, daily average air pressure, and daily evaporation, three time-varying driving factors) of each grid node are input into the principal component analysis algorithm for linear dimensionality reduction. The principal component analysis algorithm solves for the eigenvalues ​​and eigenvectors of the covariance matrix, projects the original high-dimensional feature space onto the orthogonal direction with the largest variance, and extracts the top few principal components with a cumulative variance contribution rate of 0.85 as eigenvectors.

[0020] It should be noted that a cumulative variance contribution rate threshold of 0.85 can effectively filter out low-variance components introduced by local measurement errors and random noise while retaining more than 85% of the variation information of the original data. Comparative verification shows that when the threshold is increased from 0.75 to 0.85, the number of retained principal components increases from 3 to 5, and the root mean square error of out-of-bag prediction of the random forest model decreases from 0.28m to 0.21m, with an accuracy improvement of about 25%. When the threshold is further increased from 0.85 to 0.95, the number of principal components increases to 8, but the root mean square error only decreases from 0.21m to 0.20m, with an accuracy improvement of less than 5%.

[0021] After obtaining the feature vectors, they are subjected to min-max normalization, which linearly maps each feature value in a dimension to the interval between 0 and 1 according to the maximum and minimum values ​​in that dimension, thus obtaining the input feature set. It should be noted that the maximum and minimum values ​​used for normalization are obtained statistically from the training set and are stored in a fixed manner. The normalization parameters stored in the training stage are directly used for mapping in the subsequent inference stage to ensure that the feature transformations in training and inference are consistent.

[0022] For illustrative purposes, the following normalization example only shows two principal components after PCA analysis. The principal components output by the principal component analysis are the projection values ​​of each dimension of the original geological difference data onto the eigenvector direction. The physical dimensions of the original parameters are no longer retained in the principal components, but their numerical forms can be used directly for ease of subsequent normalization. In this node, the value of the first principal component is 0.8, and the value of the second principal component is 0.03. The value range of the first principal component across the entire region is 0.1 to 1.5, and the value range of the second principal component across the entire region is 0.01 to 0.05. After normalization, the first principal component is mapped to (0.8-0.1) / (1.5-0.1)=0.50, and the second principal component is mapped to (0.03-0.01) / (0.05-0.01)=0.50, both normalized to the same numerical scale.

[0023] The normalized input feature set is used as the independent variable of the model. Historical groundwater level monitoring data from no less than 20 monitoring wells in the target area over the past 36 months (covering 3 complete hydrological years and including the complete range of water level changes during the dry and wet seasons) are used as training labels (the total number of samples is no less than 5,000 (this sample size meets the empirical requirement of at least 500 samples per dimension in a 10-dimensional feature space for the random forest model), and the samples are randomly divided into training and validation sets at a ratio of 80% and 20%, respectively). A random forest regression model is constructed.

[0024] It should be noted that in this embodiment, the number of decision trees in the random forest algorithm is set to 100. During the training process of each decision tree, approximately 63% of the samples are randomly selected with replacement from the original training set using a bootstrap sampling method to form the training subset of the tree. The remaining approximately 37% of the samples are used as out-of-bag data to evaluate the prediction accuracy of the tree. Simultaneously, at each split node of each decision tree, the square root of the total number of features is randomly selected from all geological features and rounded down as candidate split features. In this embodiment, the input features are 10-dimensional, and the number of candidate split features is floor(√10) = 3. Finally, the arithmetic mean of the prediction results of all decision trees is used as the final prediction output of the model.

[0025] After model training is complete, the input feature set of all grid nodes in the target area is input into the trained random forest model. The model outputs the groundwater level change sequence of each grid node over the next 30 days (this prediction window is aligned with the 30-day time step of the diffusion simulation in step S13, and covers at least two complete inherent cycles of water level fluctuations to ensure the statistical reliability of FFT frequency extraction). This water level change sequence is organized in the form of discrete nodes with a time step of 1 hour (the 1-hour step ensures at least 12 sampling points per cycle of the 12.5-hour fluctuation cycle), totaling 720 time steps over 30 days. Each time node records the water level elevation value at the corresponding grid location.

[0026] For example, the input feature set of 500 grid nodes in the target area is input into a random forest model, and the model outputs the water level change sequence of each node over the next 30 days. Taking the grid node where a monitoring well is located as an example, the model predicts that the water level at this node will continue to drop at a rate of 0.2m per day, and a water level elevation value is recorded every 1 hour. The initial water level elevation on day 1 is 102.3m, and the water level elevation on day 30 is 96.3m, with a cumulative drop of 6.0m over 30 days.

[0027] The water level change sequence output by each grid node is input into the Fast Fourier Transform algorithm to convert the discrete time-domain water level signal into a continuous frequency-domain signal. The frequency component with the largest amplitude in the frequency domain is extracted as the main frequency of the water level fluctuation, and the amplitude corresponding to the frequency component is recorded as the maximum amplitude.

[0028] Specifically, a discrete water level sequence of length 720 is input into a fast Fourier transform to obtain the physical amplitude corresponding to each frequency component. The frequency component with the largest amplitude is selected as the dominant frequency of water level fluctuation, and its corresponding amplitude is taken as the maximum amplitude.

[0029] The extracted dominant frequency and maximum amplitude are used as the attribute fluctuation characteristics of the grid node. These attribute fluctuation characteristics comprehensively reflect the hydraulic conductivity and water storage and release characteristics of the aquifer. A higher dominant frequency indicates that the groundwater flow system responds faster to external disturbances (such as rainfall or artificial pumping) and the hydraulic conductivity coefficient of the aquifer is larger. A larger amplitude indicates that the water level fluctuates significantly during the recharge and discharge cycles and the effective porosity or water storage coefficient of the aquifer is higher.

[0030] For example, after processing the 30-day water level change sequence of a certain grid node using Fast Fourier Transform, its frequency domain signal shows a significant peak at 0.08h⁻¹, corresponding to an amplitude of 1.5m, while the amplitudes of the other frequency components are all below 0.3m. This indicates that the aquifer at this grid node exhibits a high-frequency, large-amplitude fluctuation characteristic in response to rainfall replenishment; that is, the fluctuation period corresponding to the dominant frequency is 12.5h, and the maximum amplitude is 1.5m, indicating that the aquifer in this area has good hydraulic conductivity and a high storage-release coefficient. It should be noted that the 12.5h periodic water level fluctuation originates from the frequency domain superposition of multiple driving factors, including diurnal pressure variation, rainfall infiltration response, and solid tides. This periodic characteristic can serve as an important indicator for identifying the aquifer's sealing properties and elastic storage coefficient.

[0031] It is worth noting that, according to the statistical analysis of the attribute fluctuation characteristics of all grid nodes in the entire region, the highest dominant frequency is 0.08h⁻¹, corresponding to a minimum fluctuation period of 12.5h. The attribute fluctuation characteristics are a spatially distributed dataset with grid nodes as the unit. Different grid nodes exhibit different attribute fluctuation characteristic values ​​due to differences in stratigraphic structure and lithological combinations.

[0032] S12, acquire monitoring data fragments of the target area, determine whether the monitoring data fragments have spatiotemporal misalignment based on the attribute fluctuation characteristics, if so, reconstruct the monitoring data fragments to obtain a first monitoring sequence, if not, extract the concentration values ​​of the monitoring data fragments as a second monitoring sequence.

[0033] Specifically, the monitoring network consists of multiple monitoring wells distributed within the target area. Each monitoring well is equipped with a water level sensor and a water quality sensor. The groundwater level elevation data and the concentration data of organic pollutants in the groundwater are continuously collected according to a preset sampling frequency (in this embodiment, the sampling frequency is set once every 15 minutes. This frequency is set to 50 times the main frequency of the attribute fluctuation characteristics extracted in step S11, which is 0.08h-1 (15-minute sampling interval corresponds to a sampling frequency of 4h-1, 4 / 0.08=50 times). This far exceeds the minimum standard of twice required by the Nyquist sampling theorem. The sufficient sampling margin ensures the signal fidelity of the subsequent time-frequency analysis).

[0034] It should be noted that the interception window for the monitoring data fragments is determined by the attribute fluctuation characteristics obtained in step S11. The water level time series data and concentration time series data collected from each monitoring well within the interception window are encapsulated into an independent monitoring data fragment. Each monitoring data fragment is accompanied by a three-dimensional spatial coordinate label and a collection timestamp label for that monitoring well.

[0035] For example, in step S11, the attribute fluctuation characteristics of a certain grid node show that the period corresponding to the dominant frequency of water level fluctuation is 12.5 hours. To ensure that the interception window covers at least one complete fluctuation cycle and leaves a margin, the data fragment interception window of the monitoring wells near this node is set to 24 hours (based on the attribute fluctuation characteristic value of the grid node closest to the monitoring well; if there are multiple equidistant grid nodes, the arithmetic mean of the characteristic values ​​is taken). Within this 24-hour window, the monitoring well collects a total of 96 water level data points and 96 concentration data points at a sampling frequency of once every 15 minutes, forming the monitoring data fragment of this node.

[0036] In one implementation, determining whether the monitoring data fragments have spatiotemporal misalignment based on the attribute fluctuation characteristics includes: The fluctuation characteristics of the attribute are analyzed periodically to obtain the fluctuation period range; Based on the fluctuation period range, the monitoring data fragments are spatiotemporally aligned to obtain a fragment splicing sequence; The offset matrix is ​​calculated for the fragmented splicing sequence to obtain the offset score; If the offset score is greater than the preset offset threshold, it is determined that the monitoring data fragment has a spatiotemporal misalignment.

[0037] The reciprocal of the dominant frequency value fmain contained in the attribute fluctuation characteristics of each grid node is taken to obtain the natural period of water level fluctuation of the aquifer at that node, T = 1 / fmain. The lower limit of the fluctuation period interval is set to the natural period T, and the upper limit is set to 2T to ensure that the intercept window covers at least one complete fluctuation period.

[0038] Specifically, the timestamp and spatial coordinate labels attached to each monitoring data fragment are extracted. For data fragments with the same or adjacent timestamps, Delaunay triangulation is performed according to the spatial planar coordinates of each monitoring well. The spatial order of the data fragments is determined by the adjacency relationship of the generated triangles. Simultaneously, the timestamps are used as the third dimension for stacking, outputting a fragment splicing sequence containing spatiotemporal topological relationships. For data fragments with incompletely consistent timestamps, a time alignment tolerance of 50% of the sampling interval (i.e., 7.5 min) is set. The time-aligned data fragments are organized into a three-dimensional array structure according to the horizontal and vertical indices of the spatial coordinates. For example, the monitoring data fragments of 20 monitoring wells in the target area within a 24-hour window are spliced ​​into a fragment splicing sequence with dimensions of 20 (number of monitoring wells) × 96 (time sampling points) × 2 (bivariate variables of water level and concentration) after spatiotemporal alignment.

[0039] Cross-correlation analysis was performed on the water level change sequences of each pair of adjacent monitoring nodes in the fragmented sequence. The cross-correlation coefficient between the two sequences at different time lags was calculated, and the lag time corresponding to the maximum cross-correlation coefficient was taken as the time offset of the node pair. After performing the above calculation on all adjacent monitoring well pairs, the offsets of each node pair were organized into an offset matrix according to spatial adjacency. The mean μ of the absolute values ​​of the offsets of each node pair in the offset matrix was taken. offset With standard deviation σ offset Define the offset fraction p = μ offset +1.654×σ offset The unit is min, 1.654 × σ offse Corresponding to the 95th percentile.

[0040] Monitoring data from the past 12 months in which no known pollution events occurred in the target area were selected as the normal baseline dataset. The offset score of the normal data for each time period was calculated, and the 95th percentile was taken as the judgment threshold (approximately 15 minutes in this embodiment). When the offset score is greater than 15 minutes, it is determined that there is a spatiotemporal misalignment in the monitoring data fragments of the current batch.

[0041] It is worth noting that the physical causes of spatiotemporal misalignment mainly fall into two categories. The first is the presence of unexplored water-blocking faults in the aquifer, which interrupts or detours the hydraulic conduction path between adjacent monitoring wells, causing an abnormal delay in the time it takes for downstream monitoring wells to receive the water level response. The second is abrupt changes in local permeability coefficients, such as due to the presence of karst channels or fracture zones, where the hydraulic conduction velocity in one direction is much higher than in other directions, causing monitoring wells in that direction to detect water level changes before the expected time.

[0042] In one implementation, the step of reconstructing the monitoring data fragments to obtain a first monitoring sequence includes: Spatial relationship reorganization is performed on the monitoring data fragments with spatiotemporal misalignment to obtain a node topology map; The observation sequence is obtained by solving and filling in the state transition of the node topology graph. Spatial interpolation is performed on the observation sequence to obtain the reconstructed sequence; If the confidence level of the reconstructed sequence is greater than a preset confidence threshold, then the reconstructed sequence is smoothed and filtered to obtain the first monitoring sequence.

[0043] Each monitoring well within the target area is abstracted as a graph node, and the direction of hydraulic conduction between adjacent monitoring wells is abstracted as a directed edge, with the edge direction consistent with the groundwater flow direction. The input is a set of three-dimensional spatial coordinates of monitoring wells exhibiting misalignment characteristics. Using the Delaunay triangulation method based on the spatial planar coordinates of the monitoring wells (as described in the aforementioned spatiotemporal alignment process), spatial connections between nodes are established. The output is a node topology graph containing both the set of nodes and the set of edges.

[0044] For example, after processing by Delaunay triangulation, 20 monitoring wells with spatiotemporal misalignment in the target area are used to generate a node topology graph containing 20 graph nodes and approximately 45 edges. Each node includes the three-dimensional spatial coordinates of the monitoring well and the most recently acquired water level elevation value, and each edge includes the horizontal distance between adjacent monitoring wells and the estimated hydraulic conduction direction.

[0045] The water level status of each monitoring well in the node topology diagram is discretized into three categories over the past 24 hours: rising, falling, and stable. A water level change rate greater than 0.05 m / h is defined as a rising state, less than -0.05 m / h as a falling state, and between these two is defined as a stable state. This threshold is taken as 1.7 times the measurement uncertainty (0.03 m / h) corresponding to a 6-hour sampling interval of ±0.2 m water level sensor accuracy in the target area, ensuring that the state classification is not affected by random sensor measurement noise. For aquifers with extremely low permeability (e.g., K < 0.1 m / d), this fixed threshold may cause all data to be classified as stable. In this case, the threshold should be adaptively scaled by 10% of the amplitude of the attribute fluctuation characteristic of the nearest grid node to the monitoring well. If the monitoring well is equidistant from multiple grid nodes, the arithmetic mean of the amplitudes of these equidistant grid nodes is used as the scaling benchmark to ensure that the state classification corresponds one-to-one with the spatial distribution data output by S11. The discretized state sequence is then input into the Markov chain model. The state transition probability matrix is ​​calculated and stored based on historical monitoring data of the target area over the past 12 months (covering both the wet and dry seasons to ensure statistical representativeness). During the inference phase, the transition probability matrix calculated during the training phase is used directly, without recalculating it based on the current fragment window, to ensure data consistency between training and inference. Each element Pi→j in the matrix represents the probability of transitioning to state j at the next moment when the current node is in state i.

[0046] For time nodes with misalignment or missing data, the current state of their adjacent upstream nodes is extracted. The state of the missing node is determined based on the most likely transition direction in the state transition probability distribution. Then, the water level elevation to be filled is estimated based on the mean water level change rate corresponding to that state. For filling concentration data, the empirical regression relationship between water level and concentration is used (a linear regression c=p1×h+p0 is established on the historical synchronous monitoring data of the target area, where c is the concentration in mg / L and h is the water level elevation in m) to deduce the corresponding concentration estimate from the filled water level elevation value. In this embodiment, p1=2.3mg·L-1 / m and p0=15.7mg / L are obtained by fitting synchronous monitoring samples from the past 12 months, with a goodness of fit R2=0.76. This set of coefficients serves as an exemplary parameter for filling concentration in the target area. Other target areas can be recalibrated in the same way. After performing the above water level and concentration filling operations on all nodes with misalignment at each time step, a complete and intact observation sequence is output, which includes both water level elevation and organic pollutant concentration data.

[0047] The query state transition probability distribution matrix shows that in an upward state, there is a 60% probability of continuing the upward trend, a 30% probability of transitioning to a stable state, and a 10% probability of transitioning to a downward trend. The state with the highest probability of continuing the upward trend is taken as the predicted state for this node. Combined with the average water level change rate of 0.08 m / h for the node's historical upward trends, the estimated water level elevation at t=8 is the equivalent of increasing the elevation at t=7 by 0.08 m.

[0048] The water level elevation values ​​of each monitoring well node in the filled observation sequence are taken as known points, and the spatial grid coordinates of the target area are taken as points to be estimated (the grid resolution is the same as that in step S11, which is 50m×50m). The ordinary Kriging spatial interpolation algorithm is used for spatial expansion.

[0049] Specifically, firstly, using the spatial coordinates and water level elevations of known monitoring well nodes as input, a semivariogram function reflecting the spatial autocorrelation variation with distance is calculated. A weighted least squares method (using the number of point pairs within each distance interval as weights) is used to fit the empirical semivariogram function of a spherical model. The three parameters of the spherical model are the nugget value b0, the sill value b0+b, and the range a. The nugget value b0 represents the semivariogram value when the distance approaches zero, reflecting measurement error and microscale variation; the sill value represents the upper limit of the semivariogram function when it reaches stability; and the range a represents the maximum influence distance of spatial autocorrelation. Based on the fitted semivariogram function model, the weight coefficients of each known point for the grid node to be estimated are calculated, and the water level elevation estimate of the grid node to be estimated is calculated using a weighted linear combination. After performing the above calculations point-by-point for all blind zone grid nodes, a reconstructed sequence covering the entire target area is output.

[0050] The estimated variance reflects the degree of uncertainty in the interpolation result. A smaller value indicates more sufficient data support from known monitoring points around the node, and higher reliability of the interpolation result. The reciprocal of the estimated variance σ² is linearly mapped to an integer scalar from 0 to 100, with the confidence level m = min(100, max(0, (1 / σ² - 1 / σ²max) / (1 / σ²min - 1 / σ²max) × 100)), where σ²min and σ²max are the 5th and 95th quantiles of the reciprocals of the historical interpolation estimated variances for the target area, respectively. The mapping result is used as the confidence level of the node's value in the reconstructed sequence.

[0051] It should be noted that the preset reliability threshold is set to 85. This is based on the confidence distribution of historical batches of normal monitoring data after undergoing the same interpolation process, which has a mean of 90 and a standard deviation of 4. The threshold is set approximately 1.25 standard deviations below the mean (90 - 1.25 × 4 = 85).

[0052] When the confidence level of the reconstructed sequence is greater than 85, the sequence is input into the Savitzky-Gore filter algorithm for smoothing. When the confidence level is less than or equal to 85, the search radius of the Kriging interpolation is expanded to 1.5 times the original radius (this multiplication ensures that at least 2 to 3 nearby known monitoring wells are added to the interpolation under the condition of an average well spacing of 50m to 100m, while not exceeding the range of the semivariogram function variation 'a' to avoid introducing spatially irrelevant data), and more nearby known monitoring point data are added before re-performing the spatial interpolation. The above-mentioned radius expansion and re-interpolation operation is performed a maximum of 3 times. If the confidence level still does not meet the condition of being greater than 85 after 3 re-interpolations, the reconstructed sequence with the highest current confidence level is retained, and the uncertainty level of each node is marked. A confidence level greater than 60 and less than or equal to 85 is marked as medium uncertainty, and a confidence level less than or equal to 60 is marked as high uncertainty.

[0053] The Savitsky-Gorye filter algorithm achieves smoothing by performing low-order polynomial least-squares fitting on the data within a sliding window. In this embodiment, the filter window width is set to 7 time steps (7h), and the polynomial order is set to 3. The choice of a window width of 7h is based on the fact that this width is approximately 56% of the shortest fluctuation period of 12.5h extracted in step S11, which can effectively smooth high-frequency noise while avoiding excessive attenuation of the inherent fluctuation components of 12.5h. After the reconstructed sequence is input into the Savitsky-Gorye filter step by step, the smoothed first monitoring sequence is obtained.

[0054] For example, in the reconstructed sequence, a certain grid node experienced an abnormal water level elevation jump to 113.5m on the 15th day (the previous time was 110.2m, and the next time it recovered to 110.0m). After Savitsky-Gorye filtering, this abnormal jump was smoothed to 110.3m, which is consistent with the water level change trend before and after the time. The overall downward trend of the water level in this area (about 0.2m / d) was completely preserved.

[0055] In one implementation, extracting the concentration values ​​of the monitoring data fragments as a second monitoring sequence includes: The monitoring data fragments are checked for alignment with timestamps and spatial coordinates to obtain synchronization verification results; If the synchronization verification result meets the preset synchronization matching conditions, then the service load of the monitoring data fragment is extracted. Anomalies were removed from the service load to obtain the purity concentration value; The pure concentration values ​​are spliced ​​together in time sequence to obtain the second monitoring sequence.

[0056] The timestamp and spatial coordinate labels attached to each monitoring data fragment are directly extracted. The absolute difference between the timestamps of adjacent monitoring nodes is calculated as the time delay, and the Euclidean distance of the spatial coordinates is calculated as the spatial offset. The preset synchronization matching conditions are set as follows: the time delay does not exceed 5 minutes, and the spatial offset does not exceed 2 meters. When both indicators of adjacent monitoring nodes meet the above tolerance conditions, a synchronization verification result with a Boolean value of true is output, indicating that the batch of monitoring data fragments has good consistency in the spatiotemporal dimension and there are no undetected hidden misalignments.

[0057] It should be noted that the 5-minute time matching tolerance is set based on one-third of the sampling interval (15 minutes) of the monitoring network. The 2-meter spatial matching tolerance is set based on the GPS measurement accuracy of the wellhead coordinates (approximately 1 to 3 meters for civilian use), with the midpoint of 2 meters taken as a reasonable spatial alignment tolerance. If the synchronization verification result does not meet the preset synchronization matching conditions, it indicates that the data fragment may have a previously unidentified slight spatiotemporal misalignment. It will then be re-marked as having a spatiotemporal misalignment and transferred to the reconstruction process.

[0058] Concentration data collected by water quality sensors is extracted from monitoring data fragments as operational payloads. These payloads specifically represent the raw concentration sequences of specific organic pollutants (such as nitrates, benzene compounds, or chlorinated hydrocarbons) at the corresponding monitoring well and timestamp. For example, the raw nitrate concentration sequence obtained after extracting from a monitoring data fragment is [42, 45, 47, 43, 150, 44, 46], in mg / L.

[0059] Outlier detection is performed on the original concentration sequences in the workload using the interquartile range (IQR). Specifically, the original concentration sequences are sorted in ascending order of value, and the first quartile (Q1, 25th percentile) and the third quartile (Q3, 75th percentile) are calculated, yielding the IQR = Q3 - Q1. Anomaly detection intervals are defined as [Q1 - 1.5 × IQR, Q3 + 1.5 × IQR], with a coefficient of 1.5 being a commonly used standard for identifying outliers in statistics, independent of distribution assumptions. When a concentration value falls outside this interval, it is identified as an outlier and removed from the sequence, outputting the purified concentration value after outlier removal.

[0060] Using the timestamps corresponding to each purity concentration value as the baseline, purity concentration values ​​obtained from different monitoring frequencies and different monitoring wells are linearly stitched together along the time dimension according to the sampling time sequence. For concentration data from different monitoring wells with incompletely aligned sampling time points, linear interpolation is used to fill in the gaps on the time axis, ensuring that the stitched sequence has complete concentration values ​​on a unified time axis. The set of stitched concentration values ​​from all monitoring wells is output as the second monitoring sequence, which covers the complete organic pollutant concentration change process throughout the entire monitoring period.

[0061] S13, perform grid division and spatial interpolation on the first monitoring sequence or the second monitoring sequence to obtain the concentration distribution matrix of pollutants, and perform diffusion evolution processing on the concentration distribution matrix according to the attribute fluctuation characteristics to obtain the pollutant spread trajectory.

[0062] In one implementation, the step of performing grid division and spatial interpolation on the first monitoring sequence or the second monitoring sequence to obtain the pollutant concentration distribution matrix includes: Extract the concentration values ​​and spatial coordinates of the first monitoring sequence or the second monitoring sequence, as well as the uncertainty level label attached to each data node; The concentration values ​​and spatial coordinates are meshed to obtain spatial grid nodes; Spatial interpolation is performed on the spatial grid nodes to obtain the grid concentration values; If there are blank grids in the grid concentration value, the blank grids are locally filled to obtain the complete concentration value; The concentration distribution matrix is ​​obtained by performing matrix mapping on the complete concentration values.

[0063] Using the first or second monitoring sequence output from step S12 as input, its data structure is parsed to extract the concentration values ​​of organic pollutants and their corresponding three-dimensional spatial coordinates contained in each data record. It should be noted that the three-dimensional spatial coordinates consist of the horizontal coordinate X, the vertical coordinate Y, and the depth coordinate Z. The concentration values ​​are the measured or reconstructed concentration values ​​of the target organic pollutant (such as nitrate) at the corresponding location, in mg / L.

[0064] The spatial grid resolution is set to 10m × 10m (further refined to 10m based on the 50m grid in step S12 to improve the spatial accuracy of the diffusion simulation). The entire target monitoring area is divided into several standard square grid cells with sides of 10m on the horizontal plane. Each grid cell is uniquely identified by its row index i and column index j. In the depth direction, each planar grid cell is expanded into several layers according to the layered structure of the aquifer, with the thickness of each layer consistent with the formation thickness of the corresponding aquifer. The spatial coordinates of each monitoring node extracted in the previous step are assigned according to the grid cell they belong to, and the coordinates of the center point of the grid cell are used as the coordinates of the corresponding spatial grid node for output.

[0065] It should be noted that the 10m×10m grid resolution was set based on a comprehensive consideration of the spatial distribution density of the monitoring wells and the computational accuracy requirements of the subsequent diffusion simulation. The average spacing between the monitoring wells is approximately 50m to 100m, and the 10m resolution allows for the insertion of 4 to 9 grid nodes between adjacent monitoring wells, ensuring the smoothness and spatial continuity of the interpolation results.

[0066] Spatial interpolation is the process of extrapolating concentration values ​​from known monitoring points to grid nodes in blind areas without monitoring wells. The concentration values ​​and spatial coordinates of known monitoring nodes are used as training points, and the coordinates of all grid nodes are used as points to be estimated. A standard kriging interpolation algorithm is employed for spatial extrapolation. Specifically, the empirical semivariogram function is calculated and fitted using an exponential model, with the concentration values ​​of known monitoring points as input. The mathematical form of the exponential model is γ(r) = b0 + b[1 - exp(-r / a)], where r is the spatial distance between two points, b0 is the nugget value, b is the partial sill value, and a is the range. Based on the fitted semivariogram function model, the kriging weight coefficient wi is calculated for each grid node to be estimated, and the predicted concentration value for that grid node is calculated using a weighted linear combination. For known nodes marked as having medium uncertainty, their kriging weight is multiplied by a reduction factor of 0.7; for known nodes marked as having high uncertainty, their kriging weight is multiplied by a reduction factor of 0.3. The aforementioned reduction coefficients were determined through leave-one-out cross-validation using historical monitoring data of the target area: reduction coefficients for nodes with medium uncertainty (ranges from 0.5 to 0.9) and reduction coefficients for nodes with high uncertainty (ranges from 0.1 to 0.5) were traversed with a step size of 0.1, and the combination that minimized the root mean square error of the interpolation was selected. The validation results showed that the combinations of 0.7 and 0.3 corresponded to the lowest root mean square errors. The output is the grid concentration value covering all grid nodes.

[0067] Blank grid filling is a local secondary filling operation performed on grid cells that still have missing concentration values ​​after spatial interpolation. Under complex hydrogeological conditions, some grid cells may have missing values ​​because the surrounding known monitoring points are too sparse or they are located on the edge of the monitoring area, making it impossible for Kriging interpolation to assign effective weight coefficients to them. For such blank grids, an inverse distance weighted interpolation algorithm is used for local filling. Taking the center coordinates of the blank grid as the center, a search radius of 50m is set, and all grid nodes with known concentration values ​​within this search radius are extracted. The weight is calculated according to the reciprocal of the distance between each known node and the blank grid, with the closer the known node, the greater the weight. The weighted average is used as the filled concentration value for the blank grid. For boundary blank grids located at the edge of the monitoring area and without any known nodes within a 50m search radius, the search radius is successively expanded to 2 times, 3 times the original radius, until at least one known node is found. If no known node is found after 3 times the radius, the concentration value of the nearest effective grid node to the boundary blank grid is directly filled and marked as high uncertainty. After performing the above filling operation on all blank grids one by one, the output is the complete concentration value with all grid nodes having valid concentration values.

[0068] It should be noted that the search radius of 50m is set to match the average spacing of monitoring wells (50m to 100m), ensuring that valid data from at least 1 to 2 known monitoring wells can be covered within the search range of each blank grid. In the inverse distance weighted interpolation, the power of the distance is set to 2, meaning the weight is inversely proportional to the square of the distance. This setting allows the influence of closer known points to decay rapidly, avoiding excessive smoothing effects from distant known points on the filling results.

[0069] Matrix mapping is the operation of transforming spatially gridded complete concentration values ​​into a structured mathematical expression. Following the order of row index i (corresponding to the vertical coordinate Y) and column index j (corresponding to the horizontal coordinate X), the complete concentration values ​​of each grid node are sequentially mapped into a multidimensional numerical matrix. The value of the element ci,j in the i-th row and j-th column of the matrix is ​​the concentration value of the corresponding grid cell. If multiple depth layers are involved, the matrix is ​​expanded into a three-dimensional tensor ci,j,l, where l is the depth layer index. This matrix is ​​the pollutant concentration distribution matrix, where each element directly corresponds to the organic pollutant concentration value at a specific spatial location in the target area, providing a spatialized initial concentration field for the subsequent construction of pollution spread trajectories.

[0070] In one implementation, the step of performing diffusion evolution processing on the concentration distribution matrix based on the attribute fluctuation characteristics to obtain the pollutant spread trajectory includes: Based on the attribute fluctuation characteristics, the concentration distribution matrix is ​​corrected for terrain resistance to obtain a corrected concentration matrix; The diffusion parameters are solved by solving the modified concentration matrix to obtain the diffusion evolution parameters; The diffusion evolution parameters are used to calculate the attenuation of the modified concentration matrix to obtain the diffusion boundary profile. The spread boundary contours are overlaid in time to obtain the spread trajectory of the pollutants.

[0071] The topographic resistance correction is a step that incorporates the hindering effect of geological media on pollutant migration into the concentration distribution matrix. The attribute fluctuation characteristics obtained in step S11 not only reflect the temporal variation characteristics of groundwater level but also indirectly characterize the permeability of the aquifer medium. The amplitude parameters in the attribute fluctuation characteristics of each grid node are input into the Darcy velocity resistance model to calculate the topographic resistance value at that node. The Darcy velocity resistance model is based on the basic form of Darcy's law: q = K × J, v = q / φ, where q is the Darcy flux (specific discharge, unit m / d), K is the permeability coefficient (m / d), J is the hydraulic gradient (dimensionless), and φ is the effective porosity. The Darcy seepage velocity v = q / φ is used for solute transport calculation. A mapping relationship is established between the amplitude A in the attribute fluctuation characteristics and the permeability coefficient K: K = α × Aβ, where A is the water level fluctuation amplitude extracted in step S11, unit m; α is a dimensional calibration coefficient with dimensions (m / d) / mβ, and β is a dimensionless exponent.

[0072] The physical basis of this mapping relationship is that in unconfined aquifers, when the external recharge is similar and the storage coefficients are not significantly different, the water level response to recharge in areas with higher permeability is characterized by rapid rise and rapid discharge, with the amplitude depending on the ratio of K to the storage coefficient. This power-law positive correlation is applicable to the unconfined aquifer conditions in the target area of ​​this embodiment (the storage coefficient varies little between 0.15 and 0.30). α and β are empirical coefficients calibrated using the least squares method based on historical pumping test data of the target area. The calibration process involves selecting at least 5 pumping test wells from different hydrogeological units within the target area, recording the measured permeability Kmeas of each well and the attribute fluctuation characteristic amplitude A of the grid node where the well is located, and using the minimization of Σ(Kmeas-α×Aβ)2 as the objective function. The Levenberg-Marquardt algorithm is used to perform nonlinear fitting on α and β, and the goodness of fit R2 should be no less than 0.85 to be acceptable. In this embodiment, α=3.2, β=0.8, and the goodness of fit R²=0.91 were obtained after calibration. In other types of aquifers (such as deep confined aquifers or karst conduit flow systems), the form of this relationship may be different and needs to be recalibrated.

[0073] The concentration gradient vectors of each grid node in the corrected concentration matrix in the X, Y, and Z directions are extracted and combined with the hydraulic gradient data calculated from the water level elevation values ​​of each monitoring well in step S12, and input into the convection-diffusion equation solver. The hydraulic gradient J is the difference in water head Δh between adjacent monitoring wells divided by the distance between wells ΔL. If adjacent monitoring wells are located at different depths, they are first grouped according to the same aquifer and then the hydraulic gradient within the corresponding layer is calculated. The convection-diffusion equation describes the coupled process of convective migration and molecular diffusion-mechanical dispersion of pollutants in groundwater. Its one-dimensional form is ∂c / ∂t=D×∂2c / ∂x2-v×∂c / ∂x-λ×c, where c is the concentration, D is the hydrodynamic dispersion coefficient, q is the Darcy flux, and λ is the degradation rate constant. The value of λ was determined through an indoor tracer test on undisturbed soil columns in the target aquifer. The test temperature was controlled at 15±2℃ to represent the annual average groundwater temperature of the target aquifer (derived from 12 consecutive months of measured water temperature data from three deep monitoring wells in the target area). The test period was no less than 60 days to fully capture the initial decay stage of the first-order degradation kinetics. Alternatively, the value could be selected by referring to literature values ​​of the first-order degradation rate constant of similar organic pollutants in groundwater. In this embodiment, the value of λ for the target pollutant nitrate at a groundwater temperature of 15℃ was 0.002d⁻¹. This value was derived from 36 consecutive months of historical water quality monitoring data in the target area and calibrated using exponential decay fitting (fitted sample size n=180, R²=0.88). If local calibration data is lacking in the target area, the nitrate degradation rate constant range of 0.001d⁻¹ to 0.005d⁻¹ under similar hydrogeological conditions in the literature can be used as a reference. The solution incorporates empirical hydrogeological constraints, namely, lateral dispersion αT = 0.1 × αL and vertical dispersion αV = 0.01 × αL (this ratio is derived from Gelhar et al., Water Resources Research, 1992, 28(7): 1955-1974, a review of dispersion data from 59 sites worldwide). This reduces the dispersion parameters in the three directions to a single independent parameter, the longitudinal dispersion αL, thus transforming the underdetermined inverse problem into a well-posed problem. The longitudinal dispersion αL is determined by inverting the historical concentration field: using the historical concentration distribution matrix obtained in step S13 as the observed value and the forward-modeled concentration field of the three-dimensional convection-diffusion equation as the simulated value, the objective function E = Σ(cobs,i-csim,i)2 is constructed. The empirical value of 5m from the literature is used as the initial value of αL, and the Levenburg-Marquardt algorithm is used for iterative solution. The solution is stopped when the objective function drops below 10⁻⁴ or the number of iterations reaches 200. In this embodiment, the average value of αL obtained by inversion is 5.2m over the entire region, and the value ranges from 2.1m to 12.8m. Therefore, the average value of αT is 0.52m and the average value of αV is 0.052m. The hydrodynamic dispersion coefficient is calculated according to D=αL×v+D0, where v is the Darcy seepage velocity and D0 is the molecular diffusion coefficient.By solving the equations using finite difference discrete solutions at each grid node, the diffusion evolution parameters, including longitudinal dispersion αL, lateral dispersion αT, and vertical dispersion αV, are output.

[0074] The anisotropic dispersion and groundwater velocity data from the diffusion evolution parameters are input into the numerical solution framework of the convection-diffusion equation. The numerical solution employs an alternating directional implicit scheme, with each time step divided into two half-steps. The first step uses implicit differencing in the x-direction and explicit differencing in the y-direction; the second step uses implicit differencing in the y-direction and explicit differencing in the x-direction. Spatial discretization in both directions uses a central difference scheme. The degradation term λ×c is implicitly incorporated into the source term in each half-step to ensure numerical stability. The time step Δt (where q is the Darcy flux) is determined by the stability constraint C = q×Δt / Δx ≤ 0.5. In this embodiment, Δx = 10m (grid resolution) and Darcy flux q = K × J = 4.4 × 0.005 = 0.022m / d are taken, and Δt ≤ 227 days are calculated. In order to keep in line with the 30-day prediction window in step S11 and the subsequent propagation boundary output cycle, the diffusion evolution time step is set to 30 days to meet the above stability conditions.

[0075] Starting from the initial moment, the concentration decay value of each grid node is calculated at each time step, and the contour lines where the concentration value equals the safety threshold at each time step are extracted as the spread boundary profile. The safety threshold is set according to the Class III water standard limit for target organic pollutants in the "Groundwater Quality Standard" (GB / T 14848-2017) (nitrate is 20 mg / L as N). For example, the safety threshold for nitrate is set to 20 mg / L. In the concentration field at the 10th time step (corresponding to day 50), all grid nodes with a concentration of 20 mg / L are extracted and connected by a smooth curve to obtain the spread boundary profile at that moment. The area enclosed by this profile represents the range of the pollution plume where the nitrate concentration exceeds the safety threshold on day 50.

[0076] It should be noted that the attribute fluctuation characteristics (dominant frequency and amplitude) extracted in step S11 are considered as time-stable inherent properties of the aquifer within the simulation period. The physical basis for this assumption is that attribute fluctuation characteristics reflect the inherent response pattern of the aquifer medium skeleton to hydraulic disturbances (determined by pore structure and hydraulic conductivity), rather than short-term external driving conditions. When the aquifer structure does not change significantly during the simulation period (e.g., without large-scale artificial excavation or tectonic activity), the statistical characteristics extracted within the 30-day observation window can be reasonably extrapolated to longer simulation periods. The spread boundary contours of each time step are sequentially superimposed onto the same spatial coordinate system according to the time series. Cubic spline interpolation is used to smooth the transition between contour boundaries at adjacent times, eliminating the step-like discontinuities caused by discrete time steps. The final output is the pollutant spread trajectory, which is a sequence of spatial polygons ordered by time. Each polygon records the spatial extent, area, and average concentration of the pollution plume at the corresponding time, fully describing the positive migration and evolution process of organic pollutants in the groundwater of the target area.

[0077] For example, the temporal superposition of the nitrate pollution spread boundary contours on days 0, 50, 100, and 150 shows that the pollution plume continuously expanded in a southeast direction from its source. The affected area was approximately 2500 m² on day 50, approximately 4800 m² on day 100, and approximately 7200 m² on day 150. Simultaneously, the concentration at the center of the pollution plume gradually decreased from an initial 85 mg / L to approximately 52 mg / L. This spread trajectory visually reveals the direction of pollutant propagation, the rate of expansion, and the attenuation pattern in groundwater.

[0078] S14, extract the pollutant diffusion path based on the pollutant spread trajectory, and perform reverse tracing on the pollutant diffusion path to obtain the pollution source, including: The path skeleton of the pollutant spread trajectory is extracted to obtain the pollutant diffusion path; The connectivity and deflection angle of the pollutant diffusion path are calculated to determine whether there are branching nodes in the pollutant diffusion path. If the bifurcation node exists, then the source vector is constructed and the water catchment retention analysis is performed on the pollutant diffusion path with the bifurcation node to obtain the pollution source; If the bifurcation node does not exist, then flow tube extrapolation and isotope abundance comparison are performed on the pollutant diffusion path where the bifurcation node does not exist to obtain the pollution source.

[0079] The spread boundary contours at each time step, contained in the spread trajectory output in step S13, are input into the discrete curvature calculation model. The geometric curvature value is calculated point by point along the spread boundary contour at each time step. For a point on the planar curve, the discrete curvature is approximated by the rate of change of the angle formed by that point and its adjacent points. Points with significant changes in curvature value are extracted as geometric feature points of the contour. Geometric feature points on the boundary contour at adjacent time steps are matched and connected according to curvature similarity.

[0080] Based on the matched geometric feature point data, a Voronoi diagram-based central axis extraction algorithm is used to generate the path skeleton topology. The algorithm first discretizes the spread boundary contour into equally spaced boundary sampling points (the sampling interval is half of the 10m grid resolution in step S13, i.e., 5m, to fully capture local changes in boundary curvature). Using these sampling points as seed points, a Voronoi diagram is constructed. Line segments located inside the contour and belonging to the Voronoi boundaries of two different seed points are extracted as candidate central axis segments. All candidate segments are connected according to their endpoint proximity to form a continuous central axis. This central axis eliminates redundant morphological information caused by uneven local diffusion at the edges of the spread trajectory, retaining only the core channels reflecting the migration of pollutants along the main direction of groundwater flow. The nodes connected sequentially on the central axis are used as diffusion path nodes, outputting the pollutant diffusion path.

[0081] For example, after extracting the boundary contour of the nitrate spread trajectory along the central axis, a pollutant diffusion path extending from northwest to southeast is generated, with a total length of approximately 320m and containing 18 path nodes. The spacing between the nodes varies, with larger node spacing (approximately 25m) in areas where the pollution plume widens and smaller node spacing (approximately 10m) in areas with narrow flow channels. The density distribution of the path nodes visually reflects the changes in the width of the underground water flow channel.

[0082] The pollutant diffusion path is transformed into a spatial graph structure G=(V,E), where the node set V represents the path nodes on the diffusion path, and the edge set E represents the directed connections between adjacent path nodes. The out-degree dout(u) of each node in the graph structure is traversed, which is the number of directed edges originating from that node. Nodes with an out-degree greater than 1 indicate that at least two downstream branches have branched off from that node and are marked as candidate fork nodes.

[0083] For each node with an out-degree of 1, its unique outgoing edge direction is the downstream extension direction of that node. For candidate bifurcation nodes (out-degree ≥ 2), each outgoing edge direction represents a different branch direction. The incoming edge direction of the candidate bifurcation node is taken as the upstream flow direction, and the angle between each outgoing edge direction and the upstream flow direction is calculated. A preset bifurcation critical threshold of 30 degrees is set. This threshold is based on the physical characteristics of groundwater pollution plume bifurcation and literature experience. In heterogeneous aquifers, when a local low-permeability lens causes streamline deflection, a deflection angle within 30 degrees indicates continuous bending of the streamline (which can still be classified into the same flow tube bundle), while a deflection angle exceeding 30 degrees indicates that the streamline has undergone topological bifurcation (the branch flow tube bundle and the main flow tube bundle can be separated spatially). When a candidate bifurcation node has at least one outgoing edge direction with an angle greater than 30 degrees between it and the upstream flow direction, it is determined that a bifurcation node exists in the diffusion path. The spatial location of the bifurcation node is the coordinate of the candidate bifurcation node.

[0084] For example, after performing node out-degree analysis on the diffusion paths of 18 nodes, node 11 has an out-degree of 2 (having two downstream directed edges) and is marked as a candidate for bifurcation node; the out-degrees of the remaining nodes are all 1. The upstream flow direction of node 11 is 115 degrees southeast, and its two out-edge directions are 115 degrees southeast (main branch continuation) and 45 degrees northeast (lateral branch). The angle between the lateral branch and the upstream flow direction is 115°-45°=70°, which is greater than the critical threshold of 30 degrees, so node 11 is determined to be a bifurcation node. At the same time, node 8 also has an out-degree of 2, branching out a lateral branch (node ​​8→node 15→node 16→node 17). The direction of this branch is 90 degrees east, and the angle between this branch and the upstream flow direction (115 degrees southeast) is 25°, which is less than the critical threshold of 30 degrees, so it is not determined to be a bifurcation node. This lateral branch is included in the continuous curvature of the streamline.

[0085] The presence of bifurcation nodes indicates the existence of heterogeneous flow around the groundwater flow field. The physical cause is usually the obstruction of low-permeability lenses or water-blocking dikes, forcing the pollutant plume to split into two branches: one continues to migrate along the original main direction, while the other diffuses laterally along the bypass channel. If this bifurcation phenomenon is not identified and addressed during reverse tracing, the tracing path will extend backward along the incorrect branch, ultimately locating a false pollution source.

[0086] The reverse tracing strategy for bifurcation nodes is as follows: Concentration sampling data from monitoring wells surrounding the bifurcation node (all wells within a radius of approximately 100m, centered on the bifurcation node and with a radius twice the average spacing between wells) are processed using the same ordinary Kriging interpolation algorithm as described earlier. The output is a three-dimensional concentration gradient field for that node, with the gradient pointing in the direction of the fastest concentration increase, which is the approximate direction of the potential pollution source. Simultaneously, velocity monitoring data upstream and downstream of the bifurcation node are acquired, and the velocity decay is calculated, i.e., the percentage decrease in downstream velocity relative to upstream velocity. Velocity decay reflects the degree of energy dissipation at the node; a larger decay indicates a more significant impact of geological resistance on the node.

[0087] First, the magnitude of the concentration gradient and the velocity attenuation are normalized to both using a maximum-minimum method, mapping them uniformly to a dimensionless interval of 0 to 1. The direction of the concentration gradient is used as the direction of the source tracing vector, and the weighted average of the normalized magnitude of the concentration gradient and the normalized velocity attenuation is used as the magnitude of the source tracing vector, with each weight set to 0.5. This equal weighting is based on the assumption that the concentration gradient and velocity attenuation are equally important in determining the source tracing direction. Backtracking verification at three known pollution sources in the target area showed that the average source tracing distance error under the equal weighting setting was 12.5 m, which is better than the combinations of 0.3 / 0.7 (15.8 m) and 0.7 / 0.3 (14.2 m). For example, at a certain bifurcation node, the calculated concentration gradient is 0.5 mg / L / m, with a direction of 325 degrees northwest; the velocity attenuation is 0.2 m / d (upstream velocity 0.8 m / d, downstream velocity 0.6 m / d). The normalized value of the concentration gradient over the entire region is 0.65, the normalized value of the velocity decay is 0.55, and the magnitude of the source vector after weighted averaging is 0.5×0.65+0.5×0.55=0.60, with the direction being 325 degrees northwest.

[0088] Within each time step, virtual particles are moved in the reverse direction based on the groundwater velocity vector field at that moment (pre-simulated and stored as spatiotemporal grid data by the MODFLOW numerical model of groundwater flow in the target area. The permeability coefficient field of the MODFLOW model uses the spatial distribution K value obtained by calibration in step S13, the lateral boundary is determined by the far-field groundwater level in the regional hydrogeological survey report, the top plate receives rainwater infiltration recharge, the bottom plate is set as the impermeable boundary, the model parameters are calibrated by the measured water level data of 5 observation wells in the target area (the root mean square error of water level fitting is 0.18m), and the spatial resolution of the velocity vector field output by the simulation is consistent with the S13 grid at 10m), and the spatial trajectory points of the particles are recorded. It should be noted that the MODFLOW model is simulated offline once before the method execution based on the spatial distribution K value updated in step S13, generating a transient groundwater velocity vector field covering a 365-day tracing period. The time step is consistent with the diffusion evolution time step of step S13, which is 30 days. The velocity vector field is stored in NetCDF format, and each record contains a time index, grid row and column number, depth layer index, and velocity components in the three directions vx, vy, and vz. In the reverse tracing stage, the pre-stored data is read directly according to the time index and grid coordinates, and the MODFLOW model is not run repeatedly during particle tracing. Iteration stops when the particle trajectory converges to a certain spatial region or reaches the preset maximum tracing time. In this embodiment, the maximum tracing time is set to 365 days, which covers the seasonal groundwater flow field change cycle of a complete hydrological year (including the changes in flow direction and velocity during the wet and dry seasons), ensuring that the tracing process experiences all possible flow field states and outputs the reverse evolution trajectory. The catchment areas that this reverse evolution trajectory passes through in the spatial geographic information database are extracted, that is, the geographic units where the terrain is low-lying and groundwater collects.

[0089] The digital elevation model (DEM) data of the catchment area and the permeability coefficient of the aquifer medium are input into the hydrodynamic retardation model. The DEM provides information on the topographic relief of the catchment area, and the permeability coefficient reflects the resistance of the strata to water flow. The residence time of the trajectory in the catchment area is calculated by combining the travel length of the reverse evolution trajectory within the catchment area. The residence time τ = L × φ / (K × S), where L is the travel length of the reverse evolution trajectory in the catchment area (m), φ is the effective porosity, K is the average permeability coefficient of the aquifer in the catchment area (m / d), and S is the topographic slope (dimensionless) along the trajectory direction in the catchment area. This formula uses the solute transport seepage velocity v = q / φ = K × S / φ as the transport velocity, and L divided by the seepage velocity yields the residence time. For example, with a travel length L = 85m, an average permeability coefficient K = 4.4m / d, a topographic slope S = 0.005, and a retention time τ = 85 × 0.25 / (4.4 × 0.005) = 21.25 / 0.022 ≈ 966d ≈ 2.6 years, this retention time is much greater than the preset threshold of 100 days, indicating that the catchment area has the hydraulic conditions for long-term accumulation of pollutants.

[0090] A preset threshold for the residence time was set at 100 days. This threshold was set based on the fact that the natural decay half-life of common organic pollutants in groundwater is typically between 30 and 180 days. 100 days was chosen as the minimum residence time for moderately persistent pollutants to potentially accumulate in a localized retention area. It should be noted that the residence time calculated in the example above (approximately 966 days) is significantly greater than this threshold, indicating that the catchment area has sufficient hydraulic retention conditions for long-term pollutant accumulation, thus identifying this area as a potential pollution source retention zone. The centroid coordinates where the flow velocity in the catchment area approaches zero were extracted from the reverse evolution trajectory and used as the residence point coordinates. Using the stagnation point coordinates as the initial spatial value and the highest concentration around the stagnation point as the initial value of the release intensity, the Levenberg-Marquardt algorithm is used to perform nonlinear least squares iterations on the objective function E=Σ(ωi×(cobs,i-csim,i)2)+0.01×(Xs2+Ys2), where cobs,i is the measured concentration value of the i-th monitoring well, csim,i is the value of the alternating direction implicit difference forward modeling based on the three-dimensional convection-diffusion equation (using the concentration distribution matrix output from step S13 as the initial concentration field, with zero flux Neumann boundary and constant concentration Dirichlet boundary), and ωi is weighted according to the distance between the monitoring well and the stagnation point. The process terminates when the objective function decreases below 10⁻⁴ or the number of iterations reaches 200, and outputs the spatial coordinates of the pollution source, the release intensity, and the release start time under the first scenario.

[0091] When there are no bifurcation nodes in the diffusion path, it indicates that the pollutants migrate within a relatively homogeneous single aquifer channel, and their diffusion path can be approximated as a one-dimensional unidirectional flow tube model. The three-dimensional streamline data of the groundwater flow field are input into the flow tube generation algorithm. Using the centerline of the diffusion path as the axis and half the aquifer thickness at that location (approximately 0.5m in this embodiment) as the equivalent pipe diameter, a closed-bound unidirectional flow tube is generated. This flow tube defines the spatial envelope range of pollutant migration.

[0092] The permeability coefficient K of each flow segment within the unidirectional flow pipe is extracted (obtained from the Darcy velocity resistance model calibrated in step S13). The hydraulic gradient along the flow pipe axis, J = q / K, is calculated according to Darcy's law, where q is the Darcy flux. In this embodiment, J = 0.022 / 4.4 ≈ 0.005, resulting in minimal head loss per unit length and slow attenuation of water flow energy along the flow pipe; the total head loss for a 500m flow pipe is only approximately 2.5m. Since the head loss along the flow pipe in a homogeneous aquifer is nearly linearly distributed, a fixed step size is used for reverse sampling. The step size is twice the grid resolution (10m), i.e., 20m, to ensure at least one source sampling point is obtained within each grid cell, balancing sampling density and computational cost.

[0093] Starting from the downstream end of the unidirectional flow tube, the sampling point is moved upstream along the central axis of the tube, with a spatial coordinate point selected at intervals of one reverse calculation step. For example, a 500m long unidirectional flow tube is sampled in reverse steps of 20m, resulting in 25 sampling points from the end to the beginning of the tube. At each sampling point, the measured abundance values ​​of characteristic isotopes (such as nitrogen isotope 15N or carbon isotope 13C) of the target organic pollutants in the groundwater are extracted.

[0094] The isotope abundance measurements at each source tracing sampling point were compared point-by-point with the preset regional background baseline value. The regional background baseline value was determined by collecting water samples from uncontaminated background monitoring wells upstream of the target area and performing isotope mass spectrometry analysis. For example, the background baseline value of δ¹⁵N in the groundwater of the target area was 3‰. The abundance deviation value of each sampling point was obtained by subtracting the background baseline value of 3‰ from the δ¹⁵N measurements at each of the 25 sampling points, and then arranged in order from downstream to upstream of the flow tube to obtain the abundance deviation sequence.

[0095] The abundance deviation sequence was scanned point-by-point along the flow tube from downstream to upstream to identify segments that continuously and monotonically increase in the sequence. If a subsequence with a length of not less than 5 consecutive data points and a monotonically increasing deviation value existed (based on the statistics of abundance deviation sequences from 20 sets of historical uncontaminated background data in the target area, the maximum length of a continuously monotonically increasing segment caused by random fluctuations does not exceed 4, so 5 was taken as the threshold for distinguishing between the true isotope enrichment trend and random fluctuations), then the spatial tube segment corresponding to the subsequence was extracted as an anomalous mutation interval. The anomalous mutation interval characterizes the characteristic of continuous enrichment of isotope abundance in the tube segment in the reverse direction (towards the pollution source), which is a typical isotope fingerprint signal close to the pollution release source.

[0096] For example, in the abundance deviation sequence of 25 sampling points, the 10th to 18th data points show a trend of continuous monotonically increasing from 2‰ to 8‰, with 9 consecutive data points having a length greater than 5 and clear monotonicity. The 160m pipe segment (20m step size × 8 intervals) corresponding to these 9 sampling points is extracted as the abnormal mutation interval.

[0097] Using the coordinates of the upstream endpoint of the anomalous mutation interval as the initial spatial value and the highest concentration at that point as the initial value of the release intensity, the Levenburg-Marquardt algorithm is used to perform nonlinear least squares iterations on the objective function E=Σ(δobs,i-δsim,i)2+0.05×||(Xs,Ys,Zs)-(X0,Y0,Z0)||2, where δobs,i is the measured value of the isotopic abundance deviation at the i-th sampling point, δsim,i is the explicit finite difference forward modeling value based on the one-dimensional migration-diffusion-decay equation of the flow tube (with the upstream end of the flow tube as the constant concentration inflow boundary and the tube wall as the zero flux boundary), and (X0,Y0,Z0) is the geometric center of the anomalous mutation interval as the regularization anchor point. The process terminates when the objective function decreases below 10⁻⁴ or the number of iterations reaches 200, and outputs the spatial coordinates of the pollution source and the release intensity under the second scenario.

[0098] For example, source strength inversion iterative calculation was performed on the nine abundance deviation data within the above-mentioned 160m abnormal change interval. After 47 iterations, the objective function converged, and the pollution source was determined to be located at a position with a horizontal coordinate X of 1500m, a vertical coordinate Y of 2800m, and a depth coordinate Z of -35m, which is the spatial location of the pollution source under the second scenario.

[0099] In summary, this invention discloses a method for tracing the source of organic pollution in groundwater by integrating multi-source heterogeneous data. This invention extracts attribute fluctuation features from geological difference data using random forests and reconstructs spatiotemporally misaligned monitoring data accordingly. It then combines these attribute fluctuation features with diffusion evolution processing to construct a spread trajectory and uses a differentiated reverse tracing strategy based on bifurcation perception to locate the pollution source. This achieves efficient integration of multi-source heterogeneous geological data and monitoring data, improving the accuracy and reliability of tracing the source of organic pollution in groundwater.

[0100] The second embodiment of the present invention provides a groundwater organic pollution tracing system that integrates multi-source heterogeneous data, including: a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements the steps of the groundwater organic pollution tracing method that integrates multi-source heterogeneous data as described above.

[0101] It should be noted that the groundwater organic pollution tracing system that integrates multi-source heterogeneous data provided in this embodiment of the invention is used to execute all the process steps of the groundwater organic pollution tracing method that integrates multi-source heterogeneous data in the above embodiment. The working principles and beneficial effects of the two are one-to-one, so they will not be described again.

[0102] It should be noted that the device embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate, and the components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs. Furthermore, in the accompanying drawings of the device embodiments provided by this invention, the connection relationships between modules indicate that they have communication connections, which can be specifically implemented as one or more communication buses or signal lines. Those skilled in the art can understand and implement this without any creative effort.

[0103] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above descriptions are merely specific embodiments of the present invention and are not intended to limit the scope of protection of the present invention. In particular, it should be noted that any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention for those skilled in the art.

Claims

1. A method for tracing the source of organic pollution in groundwater by integrating multi-source heterogeneous data, characterized in that, include: Geological difference data of the target area is obtained, and random forest processing is applied to the geological difference data to obtain attribute fluctuation characteristics; Acquire monitoring data fragments of the target area, determine whether the monitoring data fragments have spatiotemporal misalignment based on the attribute fluctuation characteristics, if so, reconstruct the monitoring data fragments to obtain a first monitoring sequence, if not, extract the concentration values ​​of the monitoring data fragments as a second monitoring sequence; The first monitoring sequence or the second monitoring sequence is divided into grids and spatially interpolated to obtain the concentration distribution matrix of pollutants. The concentration distribution matrix is ​​then subjected to diffusion evolution processing based on the attribute fluctuation characteristics to obtain the pollutant spread trajectory. Based on the pollutant spread trajectory, the pollutant diffusion path is extracted, and the pollutant diffusion path is traced backward to obtain the pollution source.

2. The method for tracing the source of organic pollution in groundwater by fusing multi-source heterogeneous data according to claim 1, characterized in that, The random forest processing of the geological difference data to obtain attribute fluctuation characteristics includes: The geological difference data are subjected to feature extraction and normalization to obtain the input feature set; Random forest regression prediction is performed on the input feature set to obtain the water level change sequence; The amplitude and frequency of the water level change sequence are extracted to obtain the attribute fluctuation characteristics.

3. The method for tracing the source of organic pollution in groundwater by fusing multi-source heterogeneous data according to claim 1, characterized in that, The step of determining whether the monitoring data fragments have spatiotemporal misalignment based on the attribute fluctuation characteristics includes: The fluctuation characteristics of the attribute are analyzed periodically to obtain the fluctuation period range; Based on the fluctuation period range, the monitoring data fragments are spatiotemporally aligned to obtain a fragment splicing sequence; The offset matrix is ​​calculated for the fragmented splicing sequence to obtain the offset score; If the offset score is greater than the preset offset threshold, it is determined that the monitoring data fragment has a spatiotemporal misalignment.

4. The method for tracing the source of organic pollution in groundwater by fusing multi-source heterogeneous data according to claim 1, characterized in that, The process of reconstructing the monitoring data fragments to obtain the first monitoring sequence includes: Spatial relationship reorganization is performed on the monitoring data fragments with spatiotemporal misalignment to obtain a node topology map; The observation sequence is obtained by solving and filling in the state transition of the node topology graph. Spatial interpolation is performed on the observation sequence to obtain the reconstructed sequence; If the confidence level of the reconstructed sequence is greater than a preset confidence threshold, then the reconstructed sequence is smoothed and filtered to obtain the first monitoring sequence.

5. The method for tracing the source of organic pollution in groundwater by fusing multi-source heterogeneous data according to claim 1, characterized in that, The step of extracting the concentration values ​​of the monitoring data fragments as the second monitoring sequence includes: The monitoring data fragments are checked for alignment with timestamps and spatial coordinates to obtain synchronization verification results; If the synchronization verification result meets the preset synchronization matching conditions, then the service load of the monitoring data fragment is extracted. Anomalies were removed from the service load to obtain the purity concentration value; The pure concentration values ​​are spliced ​​together in time sequence to obtain the second monitoring sequence.

6. The method for tracing the source of organic pollution in groundwater by fusing multi-source heterogeneous data according to claim 1, characterized in that, The step of performing grid division and spatial interpolation on the first monitoring sequence or the second monitoring sequence to obtain the pollutant concentration distribution matrix includes: Extract the concentration values ​​and spatial coordinates of the first monitoring sequence or the second monitoring sequence, as well as the uncertainty level label attached to each data node; The concentration values ​​and spatial coordinates are meshed to obtain spatial grid nodes; Spatial interpolation is performed on the spatial grid nodes to obtain the grid concentration values; If there are blank grids in the grid concentration value, the blank grids are locally filled to obtain the complete concentration value; The concentration distribution matrix is ​​obtained by performing matrix mapping on the complete concentration values.

7. The method for tracing the source of organic pollution in groundwater by fusing multi-source heterogeneous data according to claim 1, characterized in that, The step of performing diffusion evolution processing on the concentration distribution matrix based on the attribute fluctuation characteristics to obtain the pollutant spread trajectory includes: Based on the attribute fluctuation characteristics, the concentration distribution matrix is ​​corrected for terrain resistance to obtain a corrected concentration matrix; The diffusion parameters are solved by solving the modified concentration matrix to obtain the diffusion evolution parameters; The diffusion evolution parameters are used to calculate the attenuation of the modified concentration matrix to obtain the diffusion boundary profile. The spread boundary contours are overlaid in time to obtain the spread trajectory of the pollutants.

8. The method for tracing the source of organic pollution in groundwater by fusing multi-source heterogeneous data according to claim 1, characterized in that, The step of extracting the pollutant diffusion path based on the pollutant spread trajectory and performing reverse tracing of the pollutant diffusion path to obtain the pollution source includes: The path skeleton of the pollutant spread trajectory is extracted to obtain the pollutant diffusion path; The connectivity and deflection angle of the pollutant diffusion path are calculated to determine whether there are branching nodes in the pollutant diffusion path. If the bifurcation node exists, then the source vector is constructed and the water catchment retention analysis is performed on the pollutant diffusion path with the bifurcation node to obtain the pollution source; If the bifurcation node does not exist, then flow tube extrapolation and isotope abundance comparison are performed on the pollutant diffusion path where the bifurcation node does not exist to obtain the pollution source.

9. A groundwater organic pollution tracing system integrating multi-source heterogeneous data, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the groundwater organic pollution tracing method according to any one of claims 1 to 8, which integrates multi-source heterogeneous data.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When executed by a processor, the program implements the groundwater organic pollution tracing method according to any one of claims 1 to 8, which integrates multi-source heterogeneous data.

Images