Method for identifying dominant factors of regional groundwater environmental risk based on grey correlation analysis

Abstract

The present application relates to the technical field of groundwater analysis, in particular to a regional groundwater environment risk leading factor identification method fusing grey correlation analysis, comprising the following steps: collecting pollution source emission timing identification intensity mutation inflection point, combining water level and permeability parameter to simulate pollution migration and dividing dynamic hydrogeological unit, calculating source-sink grey correlation according to emission stage segmentation and correcting risk interval stability, sequencing pollution source risk contribution to identify leading pollution factor.In the present application, the emission intensity mutation inflection point is identified and the dynamic boundary of the hydrogeological unit is circled, the pollutant diffusion range is accurately matched with the monitoring data in the time and space dimensions, the calculation interval is divided according to the emission characteristics and the differentiated geometric transformation is implemented, the data interval is compressed in the explosive period and the numerical distance is enlarged in the stable period, the correlation degree distribution dispersion degree is used for secondary segmentation and mapping correction, and the stability of risk level determination is ensured and the regional environment risk leading driving factor is accurately locked.

Description

Technical Field

[0001] This invention relates to the field of groundwater analysis technology, and in particular to a method for identifying the dominant factors of regional groundwater environmental risk that integrates grey relational analysis. Background Technology

[0002] The field of groundwater analysis technology refers to a collection of technologies related to monitoring, detection, evaluation, and identification of groundwater quality, quantity changes, and environmental risks. Core aspects include the analysis of pollutant composition and concentration in groundwater, the identification of natural conditions and human activities influencing groundwater, the determination of groundwater environmental risk levels, and a comprehensive analysis of groundwater environmental safety at a regional scale. This is achieved by constructing an evaluation index system, collecting multi-source monitoring data, and conducting quantitative analysis to support groundwater pollution prevention and environmental management decisions. Among these, the traditional method for identifying dominant factors of regional groundwater environmental risk using grey relational analysis involves identifying the main influencing factors in regional groundwater environmental risk research. This method first lists indicators such as pollution source emission conditions, hydrogeological parameters, groundwater extraction data, and land use or vegetation status based on experience or standards. Then, the data for each indicator are normalized or dimensionless. Finally, by calculating the grey relational coefficient and correlation degree between the indicator sequence and the groundwater environmental risk assessment result sequence, the factors are ranked to determine the dominant influencing factors that need to be focused on in groundwater environmental risk analysis.

[0003] Existing traditional grey relational analysis relies on static indicator systems and global dimensionless processing, ignoring the spatiotemporal lag characteristics of pollutant transport in groundwater flow fields. It is difficult to distinguish the differences in the impact of stable discharge periods and sudden surge periods on water quality. The uniform geometric scaling processing mode masks the specific response characteristics of high-intensity discharge intervals, resulting in reduced sensitivity in identifying sudden pollution sources. Furthermore, the lack of dynamic spatial boundary adjustment based on solute transport fronts causes a misalignment between source strength and receptor response in the spatiotemporal dimensions, leading to deviations in relational measurement calculations and weakening the accuracy of identifying the dominant driving factors of environmental risks. Summary of the Invention

[0004] To address the technical problems existing in the prior art, embodiments of the present invention provide a method for identifying the dominant factors of regional groundwater environmental risk by incorporating grey relational analysis, including the following steps: To achieve the above objectives, the present invention adopts the following technical solution: a method for identifying the dominant factors of regional groundwater environmental risk by incorporating grey relational analysis, comprising the following steps: S1: Collect instantaneous concentration values ​​of ammonia nitrogen and heavy metals and fluid flow parameters in the sewage pipes of the chemical industrial park, identify the abrupt change in emission behavior from gradual fluctuation to surge, divide the emission behavior into stable period and surge period evolution period, and generate time-varying inflection point sequence of pollution source emission intensity. S2: Call the time-varying inflection point sequence of the pollution source emission intensity, obtain the water level elevation of the monitoring well group and the permeability data of the aquifer porous medium, simulate the position of the transport front of solute in the groundwater flow field, and establish a dynamic boundary partition set of hydrogeological units; S3: Based on the dynamic boundary partition set of the hydrogeological unit, the evaluation period is divided into multiple calculation intervals according to the time-varying inflection point sequence of pollution source emission intensity. In the interval of rising emission flux, the geometric distance of the data sequence is processed by compression transformation, and in the interval of stable fluctuation, the numerical distance is processed by amplification transformation to obtain a segmented source-sink gray relational metric matrix. S4: Using the segmented source-sink grey relational metric matrix, calculate the dispersion of the correlation metric values ​​within the water quality level interval, and perform secondary segmentation on intervals with dispersion exceeding a set threshold to generate a set of stability correction indicators for risk level intervals. S5: Using the aforementioned risk level interval stability correction index set, the risk contribution of pollution sources in the differentiated hydrogeological units is ranked, the pollution source objects that occupy the first ranking in the multiple risk zones are screened, the corresponding characteristic pollutant types and emission sources are identified, and a list of regional groundwater environmental risk dominant driving factors is obtained. As a further aspect of the present invention, the time-varying inflection point sequence of pollution source emission intensity includes a timestamp of the start of emission rate surge, a fluctuation state transition identifier, and a time-series segment index; the dynamic boundary partition set of hydrogeological units includes pollutant diffusion envelope coordinates, corrected unit boundary geometric parameters, and aquifer permeable medium partition identifier; the segmented source-sink grey relational metric matrix includes interval weighted relational values, geometric similarity coefficients, and source-sink response distance matrix elements; the risk level interval stability correction index set includes a dispersion correction factor, a secondary segmentation weight coefficient, and a stability verification score; and the list of dominant driving factors for regional groundwater environmental risk includes the name of priority control pollution sources, key hazardous pollutant categories, and a risk contribution ranking index.

[0005] As a further aspect of the present invention, the specific steps of S1 are as follows: S101: Collect instantaneous concentration values ​​of ammonia nitrogen and heavy metals and fluid flow parameters in the sewage pipes of the chemical industrial park. Perform a product operation on the concentration values ​​and flow parameters at the same sampling time to obtain the emission flux value per unit time. Perform time-series alignment and indexing of the continuously collected emission flux values ​​according to the time axis. Remove redundant noise and abnormal missing values ​​in the data transmission process to construct a multi-dimensional pollutant emission rate time series dataset. S102: Call the multidimensional pollutant emission rate time series dataset, extract the emission flux values ​​of adjacent time nodes and perform first-order difference operation, calculate the ratio of the variation difference of flux values ​​at adjacent time points to the time interval to obtain the instantaneous change slope of emission flux, use a sliding window to perform weighted smoothing processing on continuous slope values, fit the overall evolution trend of emission flux in the time dimension, and generate emission flux change slope and evolution trend vector. S103: Based on the emission flux change slope and evolution trend vector, set the emission rate fluctuation judgment benchmark threshold, compare the instantaneous change slope with the fluctuation judgment benchmark threshold, filter the data segment whose slope value continuously exceeds the threshold and whose evolution trend shows a positive growth trend, lock the critical change moment node when the emission behavior changes from a smooth fluctuation state to a numerical climbing state, and generate an emission behavior state transition moment index table. S104: Based on the emission behavior state transition time index table, the recorded mutation nodes are segmented into segments along the entire time axis. The time interval before the mutation node is defined as the stable evolution period, and the time interval after the mutation node is defined as the surge evolution period. The time limit coordinates of the emission characteristic change are marked. The start time, end time and corresponding emission state identifier of the period are combined and mapped to establish a time-varying inflection point sequence of pollution source emission intensity.

[0006] As a further aspect of the present invention, the specific steps of S2 are as follows: S201: Call the time-varying inflection point sequence of the pollution source emission intensity, collect the real-time water level elevation value of the monitoring well group and the permeability parameters of the aquifer porous medium, calculate the hydraulic gradient vector between the monitoring wells, and analyze the flow velocity and direction vector of groundwater in the medium in combination with porosity data to generate the groundwater flow field velocity vector matrix. S202: Based on the groundwater flow field velocity vector matrix, the time step for solute transport calculation is set according to the time periods of the steady period and the surge period. The numerical difference iteration operation of the convection dispersion equation is performed to track the trajectory of pollutant concentration contour lines in the aquifer space, identify the critical position of the front end of the pollutant diffusion range, extract the spatial coordinate points of the transport front at the differentiated time, and obtain the coordinate sequence of the solute transport front position. S203: For the coordinate sequence of the solute transport front position, introduce the lithological boundary and geological structure distribution data obtained from hydrogeological exploration, calculate the spatial position deviation between the theoretical simulation boundary and the real-time geological water-blocking boundary, perform clipping and correction operations on the boundary coordinates of the area restricted by the geological structure, reconstruct the pollutant diffusion spatial outline that conforms to the hydrogeological conditions, divide the independent calculation areas in the differentiated time period, and establish a dynamic boundary partition set of hydrogeological units.

[0007] As a further aspect of the present invention, the specific steps of S3 are as follows: S301: Based on the dynamic boundary partition set of the hydrogeological unit, the location of groundwater quality monitoring wells within the dynamic boundary range is locked and the water quality concentration monitoring values ​​within the corresponding time window are extracted. The water quality concentration values ​​and pollutant emission rate data are aligned with the time axis. The continuous time series is cut into independent discrete data segments using abrupt change time nodes to construct a source-sink spatiotemporal coupled segmented dataset. S302: Call the source-sink spatiotemporal coupling segmented dataset, reduce the geometric longitudinal spacing of the data sequences marked as the surge evolution period, expand the numerical spacing of the data sequences marked as the stationary evolution period, reconstruct the geometric spatial distribution of the data curves in the coordinate system within the interval, and generate a nonlinear geometric scaling transformation sequence. S303: For the nonlinear geometric scaling transformation sequence, calculate the absolute difference distance between the pollution source emission curve and the groundwater quality response curve in the transformation space point by point, calculate the geometric shape correlation coefficient by combining the global peak difference and trough difference parameters of the data sequence, quantify the consistency of the changing trends of the two curves in the differentiated evolution period and spatial unit, and establish a segmented source-sink gray relational measurement matrix.

[0008] As a further aspect of the present invention, the process of using abrupt change time nodes to cut a continuous time series into independent discrete data segments specifically involves: calculating the first-order difference rate of change of the water quality concentration monitoring values ​​along the time axis, and comparing the absolute value of the first-order difference rate of change with a preset concentration fluctuation slope threshold. The time point when the absolute value of the first-order differential rate of change is greater than the concentration fluctuation slope threshold is marked as the mutation time node, and the time interval between two adjacent mutation time nodes is extracted as the independent discrete data segment. The process of calculating the geometric shape correlation coefficient by combining the global peak difference and trough difference parameters of the data sequence is as follows: traversing the nonlinear geometric scaling transformation sequence, selecting the peak value in the absolute difference distance corresponding to the calculation point as the global peak difference, and selecting the trough value as the trough difference. Set a resolution coefficient, calculate the product of the resolution coefficient and the global peak difference, and add the product to the valley difference to obtain the correlation calculation numerator; The absolute difference distance of the real-time calculated points is added to the product to obtain the correlation calculation denominator. The ratio of the correlation calculation numerator to the correlation calculation denominator is calculated, and the obtained ratio is determined as the geometric correlation coefficient.

[0009] As a further aspect of the present invention, the specific steps of S4 are as follows: S401: Call the segmented source-sink grey relational measurement matrix, classify the correlation degree values ​​into the corresponding level intervals according to the preset water quality classification standards, calculate the standard deviation and coefficient of variation of the correlation values ​​within the intervals to quantify the dispersion of the data distribution, compare the calculated dispersion values ​​with the stability judgment threshold, filter out the unstable interval ranges where the distribution dispersion continuously exceeds the threshold, and generate a set of abnormal intervals for water quality correlation dispersion. S402: Based on the set of abnormal water quality dispersion intervals, perform granular decomposition operation on the wide intervals with excessive dispersion, discretize the original unstable intervals into several continuous and evenly distributed narrow sub-intervals, and re-retrieve the corresponding pollution source emission intensity data and water quality monitoring indicators for the newly divided sub-interval time periods to obtain the source-sink correspondence subset. S403: For the source-sink correspondence subset, recalculate the source-sink correlation response values ​​under the subdivided spatiotemporal dimensions, introduce a dispersion weighted correction coefficient to normalize and adjust the correlation value deviation under the differentiated pollution level, eliminate the distortion of evaluation values ​​caused by data fluctuations within the interval, aggregate the corrected correlation values ​​and stability score parameters of the sub-interval, and establish a set of stability correction indexes for risk level intervals.

[0010] As a further aspect of the present invention, the process of quantifying the dispersion of data distribution by the standard deviation and coefficient of variation of the correlation values ​​within the calculation interval is specifically as follows: calculate the arithmetic mean of the correlation values ​​falling within the corresponding level interval, calculate the root mean square deviation of the correlation values ​​based on the arithmetic mean, and perform a square root operation to obtain the standard deviation. Calculate the ratio of the standard deviation to the arithmetic mean, determine the obtained ratio as the coefficient of variation, and use the coefficient of variation as a numerical basis for quantifying the dispersion of the data distribution; The process of performing granular decomposition on the wide range with excessive dispersion is as follows: extract the total time span of the unstable range, divide the total time span by the preset minimum resampling time step to obtain the number of segments, and cut the unstable range along the time axis at equal intervals according to the number of segments to obtain the several continuous and evenly distributed narrow sub-ranges. The process of introducing a dispersion weighted correction coefficient to normalize and adjust the deviation of the associated numerical values ​​under the differentiated pollution levels is as follows: calculate the reciprocal of the coefficient of variation and set the reciprocal as the dispersion weighted correction coefficient; calculate the product of the source-sink associated response values ​​under the subdivided spatiotemporal dimension and the dispersion weighted correction coefficient; and determine the calculated product result as the corrected associated numerical value of the sub-interval.

[0011] As a further aspect of the present invention, the specific steps of S5 are as follows: S501: Using the aforementioned risk level interval stability correction index set, extract the correction correlation values ​​within the hydrogeological unit, combine the unit sensitivity parameters to calculate the weighted contribution value of a single pollution source to the cumulative impact on the groundwater environment, sort the pollution source objects within the unit in descending order according to the value, identify the risk primary and secondary hierarchical relationship under differentiated spatial partitioning, and generate a multi-unit pollution source risk contribution ranking matrix. S502: For the multi-unit pollution source risk contribution ranking matrix, traverse the first and second ranking results of the hydrogeological unit partition, count the recurrence frequency of the target pollution source object occupying the first place of risk contribution in the differentiated partition, set multiple risk-dominant judgment logic, screen the key pollution source object IDs whose recurrence frequency meets the dominant judgment requirements, lock the source target with cross-regional risk characteristics, and obtain the cross-partition first-place risk-causing pollution source index set. S503: Call the cross-regional first-hand hazardous pollution source index set, retrieve the chemical component characteristic spectrum and physical emission pipeline location information corresponding to the target pollution source ID, parse the emission concentration characteristic pollutant category name and geographical coordinate parameters, and perform association mapping and structured assembly of pollution source name, characteristic pollutant type and risk contribution level to obtain a list of dominant driving factors of regional groundwater environmental risk.

[0012] Compared with the prior art, the advantages and positive effects of the present invention are as follows: In this invention, by identifying the inflection point of abrupt changes in emission intensity and delineating the dynamic boundaries of hydrogeological units, the pollutant diffusion range is precisely matched with monitoring data in the spatiotemporal dimension. The calculation interval is divided according to emission characteristics and differentiated geometric transformations are implemented. The data spacing is compressed during the surge period and amplified during the stable period, highlighting the closeness of the source-sink sequence in different evolution stages, enhancing the ability to capture nonlinear response relationships, and using the degree of dispersion of correlation distribution for secondary segmentation and mapping correction to eliminate the interference of data fluctuations on the evaluation results, ensuring the stability of risk level determination and accurately locking the dominant driving factors of regional environmental risk. Attached Figure Description

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

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

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

[0016] In embodiments of the present invention, words such as "exemplarily," "for example," etc., are used to indicate that something is an example, illustration, or description. Any embodiment or design described as "exemplary" in the present invention should not be construed as being more preferred or advantageous than other embodiments or designs. Specifically, the use of the word "exemplary" is intended to present the concept in a concrete manner. Furthermore, in embodiments of the present invention, the meaning expressed by "and / or" can be both, or either one.

[0017] In the embodiments of this invention, the terms "image" and "picture" may sometimes be used interchangeably. It should be noted that, without emphasizing the distinction between them, they convey the same meaning. Similarly, the terms "of," "corresponding (relevant)," and "corresponding" may sometimes be used interchangeably. It should be noted that, without emphasizing the distinction between them, they convey the same meaning.

[0018] In this embodiment of the invention, sometimes a subscript such as W1 may be written in a non-subscript form such as W1. When the difference is not emphasized, the meaning they express is the same.

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

[0020] Please see Figure 1 This invention provides a method for identifying the dominant factors of regional groundwater environmental risk by incorporating grey relational analysis, comprising the following steps: S1: Collect instantaneous concentration values ​​of ammonia nitrogen and heavy metals and fluid flow parameters in the sewage pipes of the chemical industrial park, monitor the slope and evolution trend of emission flux at adjacent moments, identify the abrupt moment when emission behavior changes from gentle fluctuation to increase, divide the evolution period into stable period and surge period according to emission behavior, mark the time limit of emission characteristic change, and generate a time-varying inflection point sequence of pollution source emission intensity. S2: Call the time-varying inflection point sequence of pollution source emission intensity, obtain the water level elevation of the monitoring well group and the permeability performance data of the aquifer porous medium, simulate the position of the transport front of solute in the groundwater flow field, delineate the spatial range of the aquifer covered by pollutant diffusion in different time periods, use hydrogeological exploration data to correct the dynamic offset of the unit boundary, and establish a dynamic boundary partition set of hydrogeological units. S3: Based on the dynamic boundary partition set of hydrogeological units, the evaluation period is divided into multiple calculation intervals according to the time-varying inflection point sequence of pollution source emission intensity. In the interval of rising emission flux, the geometric distance of the data sequence is processed by compression transformation, and in the interval of stable fluctuation, the numerical distance is processed by amplification transformation. The geometric similarity between the pollution source emission sequence and the groundwater quality monitoring sequence in the interval is calculated respectively, and the segmented source-sink gray relational metric matrix is ​​obtained. S4: Calculate the dispersion of correlation measurement values ​​within water quality level intervals using a segmented source-sink grey relational metric matrix. Perform secondary segmentation on intervals where the dispersion exceeds a set threshold. Remap the correspondence between pollution sources and water quality indicators for the newly segmented sub-intervals. Correct the correlation response values ​​under differentiated pollution levels and generate a set of stability correction indicators for risk level intervals. S5: Using the risk level interval stability correction index set, the risk contribution of pollution sources in the differentiated hydrogeological units is ranked, the pollution source objects that rank first in the multiple risk zones are screened, the corresponding characteristic pollutant types and emission sources are identified, and a list of regional groundwater environmental risk dominant driving factors is obtained. The time-varying inflection point sequence of pollution source emission intensity includes the start time stamp of emission rate surge, fluctuation state transition identifier, and time series segment index. The dynamic boundary partition set of hydrogeological units includes pollutant diffusion envelope coordinates, corrected unit boundary geometric parameters, and aquifer permeable medium partition identifier. The segmented source-sink grey relational metric matrix includes interval weighted relational values, geometric similarity coefficient, and source-sink response distance matrix elements. The risk level interval stability correction index set includes dispersion correction factor, secondary segmentation weight coefficient, and stability verification score. The list of dominant driving factors of regional groundwater environmental risk includes the name of priority control pollution source, key hazardous pollutant categories, and risk contribution ranking index.

[0021] Please see Figure 2 The specific steps of S1 are as follows: S101: Collect instantaneous concentration values ​​of ammonia nitrogen and heavy metals and fluid flow parameters in the sewage pipes of the chemical industrial park. Perform a product operation on the concentration values ​​and flow parameters at the same sampling time to obtain the emission flux value per unit time. Perform time-series alignment and indexing of the continuously collected emission flux values ​​according to the time axis. Remove redundant noise and abnormal missing values ​​in the data transmission process to construct a multi-dimensional pollutant emission rate time series dataset. The online monitoring sensor array deployed at key nodes of the main sewage discharge pipeline in the chemical industrial park was invoked, and the data acquisition frequency was set to once every 30 seconds (i.e., (minutes), simultaneously reading the concentration reading of the online ammonia nitrogen analyzer (using ion-selective electrode method). (unit: ,Right now ) and the flow velocity reading of the electromagnetic flowmeter (unit: For each sampling time The processor performs multiplication operations. To obtain the instantaneous emission flux, for example at time [time]. =10:00:00, ammonia nitrogen concentration was collected. Pipeline fluid flow rate Then calculate the emission flux at that moment. Then, according to the timestamp The ascending order will result in consecutive calculations. The numerical values ​​are stored in a time-series database, and Z-score normalization is performed on the data sequence. An anomaly detection threshold is set. This involves removing values ​​that deviate from the mean by more than three standard deviations, filling the gaps with linear interpolation of adjacent normal values, and mapping the cleaned flux data to a multidimensional vector space with time as the horizontal axis and emission rate as the vertical axis to construct a multidimensional pollutant emission rate time series dataset.

[0022] S102: Call the multidimensional pollutant emission rate time series dataset, extract the emission flux values ​​of adjacent time nodes, perform first-order difference operation, calculate the ratio of the change difference of flux values ​​at adjacent time points to the time interval to obtain the instantaneous change slope of emission flux, use a sliding window to perform weighted smoothing on the continuous slope values, fit the overall evolution trend of emission flux in the time dimension, and generate the emission flux change slope and evolution trend vector. The reading time span is Adjacent data points and ,set up Perform first-order difference operations for 1 minute. And calculate the slope of the instantaneous change in flux values ​​at adjacent time points. ,For example , Then the instantaneous slope A sliding window of length 5 is set, and a Gaussian weighted kernel is introduced to perform a smooth convolution operation on 5 consecutive slope values ​​within the window. The weight vector is set as follows. Assume the slope values ​​for five consecutive moments within the window are The smoothed slope value The smoothed slope sequence is realigned with the original time axis, and the derivative characteristics of emission flux in the time dimension are recorded point by point to generate a vector of emission flux change slope and evolution trend.

[0023] S103: Based on the slope of emission flux change and the evolution trend vector, set the emission rate fluctuation judgment benchmark threshold, compare the instantaneous change slope with the fluctuation judgment benchmark threshold, filter the data segment whose slope value continuously exceeds the threshold and whose evolution trend shows a positive growth trend, lock the critical change moment node when emission behavior changes from a smooth fluctuation state to a numerical climbing state, and generate an emission behavior state transition moment index table. Calculate the standard deviation of slope background noise Set a benchmark threshold for judging emission rate fluctuations. For example, the calculated standard deviation of background noise is Then set the baseline threshold as The program scans the values ​​in the slope and evolution trend vector of emission flux changes one by one. When a continuous... Each time point (set) The slope value of ) All meet Furthermore, when the evolution trend vector shows a positive value, the time period is determined to be the initial stage of emission surge. The first time point that meets the conditions is locked as the critical mutation time node. The specific timestamp of the node and the corresponding slope amplitude are recorded to generate an index table of emission behavior state transition time.

[0024] S104: Based on the emission behavior state transition time index table, the recorded mutation nodes are segmented into segments on the full cycle time axis. The time interval before the mutation node is defined as the stable evolution period, and the time interval after the mutation node is defined as the surge evolution period. The time limit coordinates of the emission characteristic change are marked. The start time, end time and corresponding emission state identifier of the period are combined and mapped to establish the time-varying inflection point sequence of pollution source emission intensity. Read the mutation time nodes in the emission behavior state transition time index table For example, if the mutation time is locked at 10:30:00, the entire time axis is divided into two parts with 10:30:00 as the boundary, and the time interval before the mutation time [00:00:00] is divided into two parts. 10:30:00] is marked as state identifier "0" (stationary period), and the time interval after the mutation time [10:30:00] is... 24:00:00 is marked as status identifier "1" (surge period), and the start time is extracted for each time period. End time and status indicators Construct a tuple Using structured data, we can complete the discretized state definition of continuous emission processes and establish a time-varying inflection point sequence of pollution source emission intensity.

[0025] Please see Figure 3 The specific steps of S2 are as follows: S201: Call the time-varying inflection point sequence of pollution source emission intensity, collect the real-time water level elevation values ​​of the monitoring well group and the permeability parameters of the aquifer porous medium, calculate the hydraulic gradient vector between the monitoring wells, and combine the porosity data to analyze the velocity and direction vectors of groundwater in the medium to generate the groundwater flow field velocity vector matrix. Activate the groundwater monitoring well cluster sensors within the corresponding area to simultaneously collect real-time water level elevations from wells A, B, and C. (Unit: m), combined with the effective porosity of aquifers in the geological exploration database. (Value) ) and permeability coefficient (Value) Construct a triangular mesh to calculate the hydraulic gradient between any two wells. Using a variant of Darcy's law Calculate the actual pore flow velocity, for example, the water level in well A. Well B water level The straight-line distance between the two wells Then the hydraulic gradient Actual pore flow velocity Simultaneously, based on the three-point water level method, the azimuth angle of the flow direction is analyzed, and the flow velocity magnitude and the flow direction angle are combined into a two-dimensional vector to generate the groundwater flow field velocity vector matrix.

[0026] S202: Based on the groundwater flow field velocity vector matrix, the time step for solute transport calculation is set according to the time periods of steady and rapid periods. Numerical difference iteration calculation of the convection dispersion equation is performed to track the trajectory of pollutant concentration contour lines in the aquifer space, identify the critical position of the front end of the pollutant diffusion range, extract the spatial coordinate points of the transport front at different times, and obtain the coordinate sequence of the solute transport front position. The durations of the plateau and surge periods are divided into several discrete time steps. (Set as) (Hourly), the groundwater flow field velocity vector matrix is ​​substituted into the convection dispersion equation, and the finite difference method is used to perform hourly iterative calculations on the aquifer grid nodes. Initial concentration boundary conditions are set, and the pollutant concentration in the grid at the next time step is calculated. The concentration value is tracked to the detection limit (e.g., The coordinates of the outermost grid nodes are used to connect the coordinate points to form a closed or semi-closed transport front. For each simulation step, the coordinates of key feature points on the front are extracted. For example, in the 10th hour of the surge, the coordinates of the forward's arrival were identified. Obtain the coordinate sequence of the solute transport front position.

[0027] S203: For the coordinate sequence of solute transport front positions, lithological boundary and geological structure distribution data obtained from hydrogeological exploration are introduced to calculate the spatial position deviation between the theoretical simulation boundary and the real-time geological water-impermeable boundary. The boundary coordinates of the area restricted by the geological structure are clipped and corrected to reconstruct the spatial outline of pollutant diffusion that conforms to hydrogeological conditions, divide the independent calculation areas in different time periods, and establish a dynamic boundary partition set of hydrogeological units. Import the regional hydrogeological vector map, which contains polygon coordinate data of impermeable boundaries (such as bedrock outcrops) and weakly permeable boundaries (such as clay lenses). Project the coordinate sequence of the solute transport front position onto the map coordinate system, perform spatial topology analysis, and check whether the front coordinate points fall inside the impermeable boundary polygons. For coordinate points that fall inside the boundary, perform normal correction according to the boundary tangent direction and force them to move to the edge of the geological boundary. Eliminate invalid calculation points that cross the boundary due to the theoretical homogeneity assumption. Delineate the effective pollution diffusion sector based on the corrected front position, define the sector as an independent hydrogeological unit for the current time period, and establish a dynamic boundary partition set for hydrogeological units.

[0028] Please see Figure 4 The specific steps of S3 are as follows: S301: Based on the dynamic boundary partition set of hydrogeological units, the location of groundwater quality monitoring wells within the dynamic boundary range is locked and the water quality concentration monitoring values ​​within the corresponding time window are extracted. The water quality concentration values ​​and pollutant emission rate data are aligned with the time axis. The continuous time series is cut into independent discrete data segments using abrupt change time nodes to construct a source-sink spatiotemporal coupled segmented dataset. Locate the ID number of the groundwater quality monitoring well in the area and extract the node at the moment of abrupt change. Preset time windows before and after (e.g.) The ammonia nitrogen concentration monitoring sequence within the hour was used, and the emission rate data from step S101 was retrieved and interpolated and aligned with a uniform time granularity of 5 minutes. As a dividing point, the aligned source (emission rate) and sink (water quality concentration) data pairs are cut into... (Data snippets from the stable period) and (Surge data fragment) Two independent matrices are used to construct a source-sink spatiotemporally coupled segmented dataset.

[0029] S302: Call the source-sink spatiotemporal coupling segmented dataset, reduce the geometric longitudinal spacing of the data sequences marked as the rapid evolution period, expand the numerical spacing of the data sequences marked as the stationary evolution period, reconstruct the geometric spatial distribution of the data curves in the coordinate system within the interval, and generate a nonlinear geometric scaling transformation sequence. right Logarithmic compression transformation performed on the time axis coordinates of the data matrix (in (in minutes from the moment of the abrupt change), reducing the vertical spacing of data points in geometric space during the surge period, targeting Perform a linear expansion transformation on the data matrix and its numerical axes. (Set expansion coefficient) This method amplifies the characteristics of small fluctuations and reconstructs the shape of two curves in a dimensionless coordinate system through nonlinear mapping. For example, it transforms the originally steep upward curve during the surge period into a gentle curve. The reconstructed coordinate point sequence is then combined to generate a nonlinear geometric scaling transformation sequence.

[0030] S303: For nonlinear geometric scaling transformation sequences, the absolute difference distance between the pollution source emission curve and the groundwater quality response curve in the transformation space is calculated point by point. The geometric shape correlation coefficient is calculated by combining the global peak difference and trough difference parameters of the data sequence, quantifying the consistency of the changing trends of the two curves in the differentiated evolution period and spatial unit, and establishing a segmented source-sink gray relational measurement matrix. Transformed emission curve sequence of pollution sources Groundwater quality response curve sequence Perform point-by-point scanning and calculate the corresponding time. absolute difference Find the maximum difference by traversing the entire sequence. minimum difference Set the resolution coefficient Substitute into the formula for calculating the grey relational coefficient For example, the absolute difference at a certain moment Global maximum difference minimum difference Then the correlation coefficient at that point The mean value of the correlation coefficients of points within the sequence is calculated to quantify the morphological similarity between the source and sink, and a segmented source-sink grey relational metric matrix is ​​established.

[0031] Please see Figure 5 The specific steps of S4 are as follows: S401: Call the segmented source-sink grey relational measurement matrix, classify the correlation values ​​into the corresponding level intervals according to the preset water quality classification standards, calculate the standard deviation and coefficient of variation of the correlation values ​​within the intervals to quantify the dispersion of the data distribution, compare the calculated dispersion values ​​with the stability judgment threshold, filter out the unstable interval ranges where the distribution dispersion continuously exceeds the threshold, and generate a set of abnormal intervals for water quality correlation dispersion. The correlation values ​​in the piecewise source-sink grey relational metric matrix are divided into different confidence intervals, and the standard deviation is calculated for the correlation set of the same water quality level. With coefficient of variation Set a stability threshold If the coefficient of variation of the correlation of a certain monitoring well within a certain time period is calculated as follows: ,because This indicates that the source-sink relationship is extremely unstable within a certain period and is greatly affected by external interference. The time period is marked as an unstable interval, the start and end times are recorded, and a set of abnormal intervals of water quality correlation dispersion is generated.

[0032] S402: Based on the set of abnormal intervals of water quality correlation dispersion, perform granular decomposition operation on the wide intervals with excessive dispersion, discretize the original unstable intervals into several continuous and evenly distributed narrow sub-intervals, and re-retrieve the corresponding pollution source emission intensity data and water quality monitoring indicators for the newly divided sub-interval time periods to obtain the source-sink correspondence subset. For time periods marked as unstable intervals (e.g., 2 hours in length), high-frequency granular decomposition is performed, dividing them into 12 narrow sub-intervals of 10 minutes each. For each sub-interval, a new query request is sent to the database to extract the original high-frequency emission data and water quality monitoring data within that small time period. No downsampling or smoothing is performed during this process, and the instantaneous fluctuation characteristics of the original data are completely preserved to construct a subset of micro-source-sink correspondences.

[0033] S403: For the source-sink correspondence subset, recalculate the source-sink correlation response values ​​under the subdivided spatiotemporal dimensions, introduce a dispersion weighted correction coefficient to normalize and adjust the correlation value deviation under the differentiated pollution level, eliminate the evaluation value distortion caused by data fluctuation within the interval, aggregate the corrected correlation values ​​and stability score parameters of the sub-interval, and establish a set of stability correction indexes for risk level intervals. Recalculate the grey relational coefficient for each micro-subinterval. And introduce a dispersion weighting correction coefficient. ,in For the coefficient of variation of the data in the sub-interval, perform a weighted correction operation. For example, the original correlation degree of a certain sub-interval coefficient of variation Then the correction factor Corrected correlation The corrected correlation value is then compared with the stability score of that sub-interval (by...). The data are aggregated to establish a set of stability correction indicators for risk level intervals.

[0034] Please see Figure 6 The specific steps of S5 are as follows: S501: Using the stability correction index set of risk level intervals, the correction correlation values ​​within the hydrogeological unit are extracted. Combined with the unit sensitivity parameters, the weighted contribution value of a single pollution source to the cumulative impact on the groundwater environment is calculated. The pollution source objects within the unit are sorted in descending order according to the value, the primary and secondary risk hierarchy under the differentiated spatial partitioning is identified, and a multi-unit pollution source risk contribution ranking matrix is ​​generated. For a specific hydrogeological unit, extract the corrected source-sink correlation values. Combined with the geological sensitivity parameters of this unit (Determined by the aquifer's permeability and antifouling performance, the value ranges from 0 to 1; for example, for gravel layers...) ), calculate the risk contribution value of a single pollution source For example, for Source_B, the degree of relevance Sensitivity Then the contribution value Potential pollution sources within the computing unit The values ​​are calculated and sorted from largest to smallest, as shown in Table 1. This identifies the main risk sources under the spatial partition and generates a multi-unit pollution source risk contribution ranking matrix.

[0035] Table 1: Example Table of Pollution Source Risk Contribution Ranking

[0036] S502: For the multi-unit pollution source risk contribution ranking matrix, traverse the first and second ranking results of the hydrogeological unit partition, count the recurrence frequency of the target pollution source object occupying the first place of risk contribution in the differentiated partition, set multiple risk-dominant judgment logic, screen the key pollution source object IDs whose recurrence frequency meets the dominant judgment requirements, lock the source target with cross-regional risk characteristics, and obtain the first risk-causing pollution source index set across partitions. Traverse the 10 different hydrogeological zones within the park and count the number of times each pollution source ranks first (Rank=1) in each zone. The multi-risk-dominated judgment logic is set as follows: if (That is, if it is the dominant risk source in at least 3 partitions), then the pollution source is determined to have a wide-area diffusion risk. For example, if Source_A ranks first in partitions 1, 3, and 5, it meets the requirement. Based on the judgment criteria, Source_A is identified as the key risk target, and the first and second risky pollution source index set across partitions is obtained.

[0037] S503: Call the first and second risky pollution source index set across the region, retrieve the chemical component feature spectrum and physical emission pipeline location information corresponding to the target pollution source ID, parse the emission concentration feature pollutant category name and geographic coordinate parameters, and perform association mapping and structured assembly of pollution source name, feature pollutant type and risk contribution level to obtain a list of dominant driving factors of regional groundwater environmental risk. Using the ID (e.g., Source_A) in the first-to-last hazardous pollution source index set across zones, retrieve the detailed file corresponding to that ID in the enterprise's environmental management database, extract characteristic pollutant fingerprint information (e.g., "aniline" and "nitrobenzene"), and the precise GPS coordinates of the discharge outlet. The pollution source name "No. 3 Synthesis Workshop Discharge Outlet", the characteristic pollutant "aniline", and the calculated highest risk contribution value "0.680" are structurally bound together to output a complete risk tracing report and obtain a list of dominant driving factors of regional groundwater environmental risk.

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

Claims

1. A method for identifying dominant factors of regional groundwater environmental risk integrating grey relational analysis, characterized in that, Includes the following steps: S1: Collect instantaneous concentration values ​​of ammonia nitrogen and heavy metals and fluid flow parameters in the sewage pipes of the chemical industrial park, identify the abrupt change in emission behavior from gradual fluctuation to surge, divide the emission behavior into stable period and surge period evolution period, and generate time-varying inflection point sequence of pollution source emission intensity. S2: Call the time-varying inflection point sequence of the pollution source emission intensity, obtain the water level elevation of the monitoring well group and the permeability data of the aquifer porous medium, simulate the position of the transport front of solute in the groundwater flow field, and establish a dynamic boundary partition set of hydrogeological units; S3: Based on the dynamic boundary partition set of the hydrogeological unit, the evaluation period is divided into multiple calculation intervals according to the time-varying inflection point sequence of pollution source emission intensity. In the interval of rising emission flux, the geometric distance of the data sequence is processed by compression transformation, and in the interval of stable fluctuation, the numerical distance is processed by amplification transformation to obtain a segmented source-sink gray relational metric matrix. S4: Using the segmented source-sink grey relational metric matrix, calculate the dispersion of the correlation metric values ​​within the water quality level interval, and perform secondary segmentation on intervals with dispersion exceeding a set threshold to generate a set of stability correction indicators for risk level intervals. S5: Using the aforementioned risk level interval stability correction index set, the risk contribution of pollution sources in differentiated hydrogeological units is ranked, the pollution source objects that occupy the first ranking in the multiple risk zones are screened, the corresponding characteristic pollutant types and emission sources are identified, and a list of regional groundwater environmental risk dominant driving factors is obtained.

2. The method for identifying dominant factors of regional groundwater environmental risk by incorporating grey relational analysis according to claim 1, characterized in that, The time-varying inflection point sequence of pollution source emission intensity includes the start time stamp of emission rate surge, fluctuation state transition identifier, and time-series segment index. The dynamic boundary partition set of hydrogeological units includes pollutant diffusion envelope coordinates, corrected unit boundary geometric parameters, and aquifer permeable medium partition identifier. The segmented source-sink grey relational metric matrix includes interval weighted relational values, geometric similarity coefficient, and source-sink response distance matrix elements. The risk level interval stability correction index set includes dispersion correction factor, secondary segmentation weight coefficient, and stability verification score. The list of dominant driving factors for regional groundwater environmental risk includes priority control pollution source names, key hazardous pollutant categories, and risk contribution ranking index.

3. The method for identifying dominant factors of regional groundwater environmental risk by incorporating grey relational analysis according to claim 1, characterized in that, The specific steps of S1 are as follows: S101: Collect instantaneous concentration values ​​of ammonia nitrogen and heavy metals and fluid flow parameters in the sewage pipes of the chemical industrial park. Perform a product operation on the concentration values ​​and flow parameters at the same sampling time to obtain the emission flux value per unit time. Perform time-series alignment and indexing of the continuously collected emission flux values ​​according to the time axis. Remove redundant noise and abnormal missing values ​​in the data transmission process to construct a multi-dimensional pollutant emission rate time series dataset. S102: Call the multidimensional pollutant emission rate time series dataset, extract the emission flux values ​​of adjacent time nodes and perform first-order difference operation, calculate the ratio of the variation difference of flux values ​​at adjacent time points to the time interval to obtain the instantaneous change slope of emission flux, use a sliding window to perform weighted smoothing processing on continuous slope values, fit the overall evolution trend of emission flux in the time dimension, and generate emission flux change slope and evolution trend vector. S103: Based on the emission flux change slope and evolution trend vector, set the emission rate fluctuation judgment benchmark threshold, compare the instantaneous change slope with the fluctuation judgment benchmark threshold, filter the data segment whose slope value continuously exceeds the threshold and whose evolution trend shows a positive growth trend, lock the critical change moment node when the emission behavior changes from a smooth fluctuation state to a numerical climbing state, and generate an emission behavior state transition moment index table. S104: Based on the emission behavior state transition time index table, the recorded mutation nodes are segmented into segments along the entire time axis. The time interval before the mutation node is defined as the stable evolution period, and the time interval after the mutation node is defined as the surge evolution period. The time limit coordinates of the emission characteristic change are marked. The start time, end time and corresponding emission state identifier of the period are combined and mapped to establish a time-varying inflection point sequence of pollution source emission intensity.

4. The method for identifying dominant factors of regional groundwater environmental risk by incorporating grey relational analysis according to claim 3, characterized in that, The specific steps of S2 are as follows: S201: Call the time-varying inflection point sequence of the pollution source emission intensity, collect the real-time water level elevation value of the monitoring well group and the permeability parameters of the aquifer porous medium, calculate the hydraulic gradient vector between the monitoring wells, and analyze the flow velocity and direction vector of groundwater in the medium in combination with porosity data to generate the groundwater flow field velocity vector matrix. S202: Based on the groundwater flow field velocity vector matrix, the time step for solute transport calculation is set according to the time periods of the steady period and the surge period. The numerical difference iteration operation of the convection dispersion equation is performed to track the trajectory of pollutant concentration contour lines in the aquifer space, identify the critical position of the front end of the pollutant diffusion range, extract the spatial coordinate points of the transport front at the differentiated time, and obtain the coordinate sequence of the solute transport front position. S203: For the coordinate sequence of the solute transport front position, introduce the lithological boundary and geological structure distribution data obtained from hydrogeological exploration, calculate the spatial position deviation between the theoretical simulation boundary and the real-time geological water-blocking boundary, perform clipping and correction operations on the boundary coordinates of the area restricted by the geological structure, reconstruct the pollutant diffusion spatial outline that conforms to the hydrogeological conditions, divide the independent calculation areas in the differentiated time period, and establish a dynamic boundary partition set of hydrogeological units.

5. The method for identifying dominant factors of regional groundwater environmental risk by incorporating grey relational analysis according to claim 4, characterized in that, The specific steps for S3 are as follows: S301: Based on the dynamic boundary partition set of the hydrogeological unit, the location of groundwater quality monitoring wells within the dynamic boundary range is locked and the water quality concentration monitoring values ​​within the corresponding time window are extracted. The water quality concentration values ​​and pollutant emission rate data are aligned with the time axis. The continuous time series is cut into independent discrete data segments using abrupt change time nodes to construct a source-sink spatiotemporal coupled segmented dataset. S302: Call the source-sink spatiotemporal coupling segmented dataset, reduce the geometric longitudinal spacing of the data sequences marked as the surge evolution period, expand the numerical spacing of the data sequences marked as the stationary evolution period, reconstruct the geometric spatial distribution of the data curves in the coordinate system within the interval, and generate a nonlinear geometric scaling transformation sequence. S303: For the nonlinear geometric scaling transformation sequence, calculate the absolute difference distance between the pollution source emission curve and the groundwater quality response curve in the transformation space point by point, calculate the geometric shape correlation coefficient by combining the global peak difference and trough difference parameters of the data sequence, quantify the consistency of the changing trends of the two curves in the differentiated evolution period and spatial unit, and establish a segmented source-sink gray relational measurement matrix.

6. The method for identifying dominant factors of regional groundwater environmental risk by incorporating grey relational analysis according to claim 5, characterized in that, The process of using abrupt change moment nodes to cut a continuous time series into independent discrete data segments specifically involves: calculating the first-order difference rate of change of the water quality concentration monitoring values ​​along the time axis, and comparing the absolute value of the first-order difference rate of change with a preset concentration fluctuation slope threshold. The time point when the absolute value of the first-order differential rate of change is greater than the concentration fluctuation slope threshold is marked as the mutation time node, and the time interval between two adjacent mutation time nodes is extracted as the independent discrete data segment. The process of calculating the geometric shape correlation coefficient by combining the global peak difference and trough difference parameters of the data sequence is as follows: traversing the nonlinear geometric scaling transformation sequence, selecting the peak value in the absolute difference distance corresponding to the calculation point as the global peak difference, and selecting the trough value as the trough difference. Set a resolution coefficient, calculate the product of the resolution coefficient and the global peak difference, and add the product to the valley difference to obtain the correlation calculation numerator; The absolute difference distance of the real-time calculated points is added to the product to obtain the correlation calculation denominator. The ratio of the correlation calculation numerator to the correlation calculation denominator is calculated, and the obtained ratio is determined as the geometric correlation coefficient.

7. The method for identifying dominant factors of regional groundwater environmental risk by incorporating grey relational analysis according to claim 5, characterized in that, The specific steps of S4 are as follows: S401: Call the segmented source-sink grey relational measurement matrix, classify the correlation degree values ​​into the corresponding level intervals according to the preset water quality classification standards, calculate the standard deviation and coefficient of variation of the correlation values ​​within the intervals to quantify the dispersion of the data distribution, compare the calculated dispersion values ​​with the stability judgment threshold, filter out the unstable interval ranges where the distribution dispersion continuously exceeds the threshold, and generate a set of abnormal intervals for water quality correlation dispersion. S402: Based on the set of abnormal water quality dispersion intervals, perform granular decomposition operation on the wide intervals with excessive dispersion, discretize the original unstable intervals into several continuous and evenly distributed narrow sub-intervals, and re-retrieve the corresponding pollution source emission intensity data and water quality monitoring indicators for the newly divided sub-interval time periods to obtain the source-sink correspondence subset. S403: For the source-sink correspondence subset, recalculate the source-sink correlation response values ​​under the subdivided spatiotemporal dimensions, introduce a dispersion weighted correction coefficient to normalize and adjust the correlation value deviation under the differentiated pollution level, eliminate the distortion of evaluation values ​​caused by data fluctuations within the interval, aggregate the corrected correlation values ​​and stability score parameters of the sub-interval, and establish a set of stability correction indexes for risk level intervals.

8. The method for identifying dominant factors of regional groundwater environmental risk by incorporating grey relational analysis according to claim 7, characterized in that, The process of quantifying the dispersion of data distribution by the standard deviation and coefficient of variation of the correlation values ​​within the calculation interval is as follows: calculate the arithmetic mean of the correlation values ​​falling within the corresponding level interval, calculate the root mean square deviation of the correlation values ​​based on the arithmetic mean, and perform a square root operation to obtain the standard deviation. Calculate the ratio of the standard deviation to the arithmetic mean, determine the obtained ratio as the coefficient of variation, and use the coefficient of variation as a numerical basis for quantifying the dispersion of the data distribution; The process of performing granular decomposition on the wide range with excessive dispersion is as follows: extract the total time span of the unstable range, divide the total time span by the preset minimum resampling time step to obtain the number of segments, and cut the unstable range along the time axis at equal intervals according to the number of segments to obtain the several continuous and evenly distributed narrow sub-ranges. The process of introducing a dispersion weighted correction coefficient to normalize and adjust the deviation of the associated numerical values ​​under the differentiated pollution levels is as follows: calculate the reciprocal of the coefficient of variation and set the reciprocal as the dispersion weighted correction coefficient; calculate the product of the source-sink associated response values ​​under the subdivided spatiotemporal dimension and the dispersion weighted correction coefficient; and determine the calculated product result as the corrected associated numerical value of the sub-interval.

9. The method for identifying dominant factors of regional groundwater environmental risk by incorporating grey relational analysis according to claim 7, characterized in that, The specific steps of S5 are as follows: S501: Using the aforementioned risk level interval stability correction index set, extract the correction correlation values ​​within the hydrogeological unit, combine the unit sensitivity parameters to calculate the weighted contribution value of a single pollution source to the cumulative impact on the groundwater environment, sort the pollution source objects within the unit in descending order according to the value, identify the risk primary and secondary hierarchical relationship under differentiated spatial partitioning, and generate a multi-unit pollution source risk contribution ranking matrix. S502: For the multi-unit pollution source risk contribution ranking matrix, traverse the first and second ranking results of the hydrogeological unit partition, count the recurrence frequency of the target pollution source object occupying the first place of risk contribution in the differentiated partition, set multiple risk-dominant judgment logic, screen the key pollution source object IDs whose recurrence frequency meets the dominant judgment requirements, lock the source target with cross-regional risk characteristics, and obtain the cross-partition first-place risk-causing pollution source index set. S503: Call the cross-regional first-hand hazardous pollution source index set, retrieve the chemical component characteristic spectrum and physical emission pipeline location information corresponding to the target pollution source ID, parse the emission concentration characteristic pollutant category name and geographical coordinate parameters, and perform association mapping and structured assembly of pollution source name, characteristic pollutant type and risk contribution level to obtain a list of dominant driving factors of regional groundwater environmental risk.

Images