Method for evaluating pollution risk of regional scale drug pollutants to underground water

Abstract

The invention discloses a pollution risk assessment method for underground water by regional scale drug pollutants, and relates to the field of environmental engineering. The method comprises the following steps: acquiring various soil characteristic parameter data and soil profile data; inputting the various soil characteristic parameter data into a solid-liquid distribution coefficient prediction model, and analyzing a nonlinear relationship among the solid-liquid distribution coefficient of the drug pollutants, the various soil characteristic parameter data and the environmental conditions through the solid-liquid distribution coefficient prediction model to obtain the solid-liquid distribution coefficient of the drug pollutants; processing the soil profile data through an unsaturated zone attenuation coefficient prediction model to obtain a drug pollutant unsaturated zone attenuation factor; obtaining an organic pollutant underground water dilution factor of the research area; and determining the risk value of the drug pollutants through the solid-liquid distribution coefficient of the drug pollutants, the unsaturated zone attenuation factor of the drug pollutants and the underground water dilution factor of the organic pollutants. According to the method, the accuracy of determining the risk value of the drug pollutants can be improved.

Description

Technical Field

[0001] The present invention relates to the field of environmental engineering, and particularly relates to a method for evaluating the pollution risk of pharmaceutical pollutants to groundwater at the regional scale. Background Art

[0002] In recent years, the pollution of new pollutants such as pharmaceuticals and personal care products (PPCPs) to the environment has attracted much attention, and the supervision of new pollutants has become increasingly strict. Among them, pharmaceutical new pollutants cover agents used to treat human and animal diseases, including various drugs such as antidepressants, non-steroidal anti-inflammatory drugs (NSAIDs), sedatives, painkillers, various anti-inflammatory drugs, antihypertensive drugs, anticancer drugs, antipsychotic drugs, antibacterial agents, oral contraceptives, lipid regulators, synthetic hormones, antibiotics, etc.

[0003] Researchers used the vadose zone solute transport model HYDRUS-1D coupled with the groundwater dilution model to construct the soil risk control value based on protecting groundwater at the plot scale. Although the researchers' method considered the migration process of pollutants, the existing research on PPCPs pollution mainly focuses on the water environment, resulting in a lack of systematicness in the PPCPs pollution control strategy and low accuracy of the obtained soil risk values. Summary of the Invention

[0004] Based on this, it is necessary to provide a method for evaluating the pollution risk of pharmaceutical pollutants to groundwater at the regional scale for the above technical problems. This method can improve the accuracy of determining the risk value of pharmaceutical pollutants.

[0005] The present invention adopts the following technical solutions: Obtain various soil characteristic parameter data and soil profile data of the research area; the soil characteristic parameter data includes solid-phase physical and chemical properties, adsorption experimental conditions, and molecular descriptors; the soil profile data includes meteorological data and the physical and chemical properties of target PPCPs; the meteorological data includes precipitation, temperature, humidity, and evaporation; the physical and chemical properties of target PPCPs include half-life, Henry's constant, diffusion coefficient, and organic carbon adsorption coefficient; Input various soil characteristic parameter data into the solid-liquid distribution coefficient prediction model, and analyze the non-linear relationship between the solid-liquid distribution coefficient of pharmaceutical pollutants and various soil characteristic parameter data and environmental conditions through the solid-liquid distribution coefficient prediction model to obtain the solid-liquid distribution coefficient of pharmaceutical pollutants; Process the soil profile data through the unsaturated zone attenuation coefficient prediction model to obtain the unsaturated zone attenuation factor of pharmaceutical pollutants; Obtain the groundwater dilution factor of organic pollutants in the research area; Determine the risk value of pharmaceutical pollutants through the solid-liquid distribution coefficient of pharmaceutical pollutants, the attenuation factor of pharmaceutical pollutants in the unsaturated zone, and the groundwater dilution factor of organic pollutants.

[0006] Preferably, the construction process of the solid-liquid distribution coefficient prediction model specifically includes: Obtain sample data of various soil characteristic parameters, and divide the sample data of various soil characteristic parameters into a training set and a test set; Train the random forest method, the extreme gradient boosting method, and the support vector regression method respectively through the training set to obtain the hyperparameter combinations of the random forest method, the extreme gradient boosting method, and the support vector regression method; Detect the prediction accuracy of the random forest method, the extreme gradient boosting method, and the support vector regression method based on the hyperparameter combinations respectively through the test set; Determine the method corresponding to the highest prediction accuracy as the solid-liquid distribution coefficient prediction model.

[0007] Preferably, process the soil profile data through the attenuation coefficient prediction model in the unsaturated zone to obtain the attenuation factor of pharmaceutical pollutants in the unsaturated zone, specifically including: Input the soil profile data into the one-dimensional hydrological model for simulation to obtain the pollutant concentration in the pore water at a specific depth of the soil; Calculate the attenuation factor of pharmaceutical pollutants in the unsaturated zone based on the maximum value of the pollutant concentration in the pore water at a specific depth of the soil.

[0008] Preferably, the calculation method of the attenuation factor of pharmaceutical pollutants in the unsaturated zone is: ; Wherein, AF is the attenuation factor of pharmaceutical pollutants in the unsaturated zone, is the pollutant concentration in the pore water of the soil surface layer, C d is the maximum value of the pollutant concentration in the pore water at a specific depth of the soil.

[0009] Preferably, the calculation method of the groundwater dilution factor of organic pollutants is: ; Wherein, DF is the groundwater dilution factor of organic pollutants, Cgw is the concentration of organic pollutants in the groundwater after dilution, C d is the maximum value of the pollutant concentration in the pore water at a specific depth of the soil.

[0010] Preferably, the calculation method of the groundwater dilution factor of organic pollutants is: ; Among them, K is the hydraulic conductivity of the aquifer, i is the hydraulic gradient, d is the depth of the mixing zone, I is the infiltration rate, L is the length of the source area parallel to the groundwater flow.

[0011] Preferably, the calculation method of the risk value of pharmaceutical pollutants is: ; Among them, is the risk value of pharmaceutical pollutants, is the groundwater concentration limit value of pharmaceuticals, K d is the solid-liquid distribution coefficient of pharmaceutical pollutants, AF is the attenuation factor of pharmaceutical pollutants in the unsaturated zone, DF is the groundwater dilution factor of organic pollutants.

[0012] The present invention provides a device for evaluating the pollution risk of pharmaceutical pollutants to groundwater at the regional scale, including: An acquisition module, configured to acquire various soil characteristic parameter data and soil profile data of the research area; the soil characteristic parameter data includes solid-phase physical and chemical properties, adsorption experimental conditions, and molecular descriptors; the soil profile data includes meteorological data and the physical and chemical properties of target PPCPs; the meteorological data includes precipitation, temperature, humidity, and evaporation; the physical and chemical properties of target PPCs include half-life, Henry's constant, diffusion coefficient, and organic carbon adsorption coefficient; A first determination module, configured to input various soil characteristic parameter data into the solid-liquid distribution coefficient prediction model, and analyze the non-linear relationship between the solid-liquid distribution coefficient of pharmaceutical pollutants and various soil characteristic parameter data and environmental conditions through the solid-liquid distribution coefficient prediction model to obtain the solid-liquid distribution coefficient of pharmaceutical pollutants; A second determination module, configured to process the soil profile data through the unsaturated zone attenuation coefficient prediction model to obtain the attenuation factor of pharmaceutical pollutants in the unsaturated zone; A third determination module, configured to acquire the groundwater dilution factor of organic pollutants in the research area; A fourth determination module, configured to determine the risk value of pharmaceutical pollutants through the solid-liquid distribution coefficient of pharmaceutical pollutants, the attenuation factor of pharmaceutical pollutants in the unsaturated zone, and the groundwater dilution factor of organic pollutants.

[0013] The present invention provides a computer-readable storage medium, and the storage medium stores a computer program, and when the computer program is executed by a processor, it implements the above-mentioned method for evaluating the pollution risk of pharmaceutical pollutants to groundwater.

[0014] The present invention provides a computer device, including a memory, a processor, and a computer program stored on the memory and executable on the processor. When the processor executes the program, the above-mentioned method for evaluating the pollution risk of pharmaceutical pollutants to groundwater at the regional scale is implemented.

[0015] The above at least one technical solution adopted by the present invention can achieve the following beneficial effects: constructing a prediction model for the solid-liquid partition coefficient of pharmaceutical pollutants, and using various soil characteristic parameter data as model input parameters to predict the solid-liquid partition coefficient of pharmaceutical pollutants, so as to obtain the solid-liquid partition coefficient of pharmaceutical pollutants; through a machine learning algorithm, based on soil profile data, meteorological data, and the properties of PPCPs, constructing a prediction model for the attenuation coefficient in the unsaturated zone, and using the prediction model for the attenuation coefficient in the unsaturated zone to make predictions, so as to obtain the attenuation factor in the unsaturated zone; obtaining the groundwater dilution factor of organic pollutants in the study area; determining the risk value of pharmaceutical pollutants through the solid-liquid partition coefficient, attenuation factor in the unsaturated zone, and groundwater dilution factor of pharmaceutical pollutants. Determining the risk value of pharmaceutical pollutants through the three parameters takes into account the influence of the three parameters on the risk value of pharmaceutical pollutants and improves the accuracy of determining the risk value of pharmaceutical pollutants. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] The drawings described herein are used to provide a further understanding of the present invention and constitute a part of the present invention. The schematic embodiments of the present invention and their descriptions are used to explain the present invention and do not constitute an improper limitation to the present invention. In the drawings: Figure 1 It is a technical route diagram of a method for evaluating the pollution risk of pharmaceutical pollutants to groundwater at the regional scale provided by the present invention; Figure 2 It is a schematic flow diagram of a method for evaluating the pollution risk of pharmaceutical pollutants to groundwater at the regional scale provided by the present invention; Figure 3 It is a method framework diagram of a prediction model for the attenuation coefficient in the unsaturated zone provided by the present invention; Figure 4 It is a schematic diagram of the study area and sampling points provided by the present invention; Figure 5 It is a computational conceptual model diagram of HYDRUS-1D provided by the present invention; Figure 6 It is a schematic diagram of the attenuation factors in the unsaturated zone of 3 PPCPs at different groundwater depths in the Yangtze River Delta region provided by the present invention; Figure 7 It is a schematic diagram of the attenuation factors in the unsaturated zone of 3 PPCPs provided by the present invention; Figure 8 It is a schematic diagram of a device for evaluating the pollution risk of pharmaceutical pollutants to groundwater provided by the present invention; Figure 9 Schematic diagram of a computer device for implementing a method for assessing the pollution risk of pharmaceutical pollutants to groundwater provided by the present invention. Detailed implementation manners

[0017] To make the objectives, technical solutions and advantages of the present invention clearer, the technical solutions of the present invention will be clearly and completely described below in conjunction with specific embodiments of the present invention and the corresponding drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without any creative work belong to the scope of protection of the present invention.

[0018] Devices such as desktop computers, servers, and laptop computers that can execute the solutions of the present invention. For the convenience of description, only a server will be used as the execution subject for illustration below.

[0019] Taking the study area as an example, the present invention uses trimethoprim, nalidixic acid, and niflumic acid as target pollutants to evaluate the soil pollution risks of three types of pharmaceutical pollutants for protecting groundwater safety in the study area. Trimethoprim (TMP) is a common nitrogen-containing antibiotic and one of the 14 high-risk drugs with a risk coefficient higher than 10 in hospital wastewater. Clinically, it is often used in combination with sulfamethoxazole to enhance the antibacterial effect. In the past thirty years, many researchers have paid great attention to the migration and fate of trimethoprim in the environment. Researchers evaluated the dissipation half-lives of 7 PPCPs in 13 soils and found that the average dissipation half-life of trimethoprim in soils with higher quality (well-developed structure and high nutrient content) was lower, inferring that trimethoprim had a greater migration potential in the soil-groundwater environment. The drug residues nalidixic acid (NA) and niflumic acid (NFA) have attracted wide attention as new pollutants in the environment. The environmental behavior of drugs, including migration and transformation pathways, is particularly important for environmental monitoring, risk assessment, and research on the impact on human health. Researchers found through research that nalidixic acid and niflumic acid can pollute the groundwater environment by infiltrating or leaching through the soil into the groundwater.

[0020] Current research on PPCPs pollution mainly focuses on the water environment. The research on the environmental risks of PPCPs in the soil - groundwater system as a whole is relatively weak, resulting in a relatively weak scientific basis for constructing methods for assessing the environmental risks of soil PPCPs polluting groundwater. At the same time, the current soil environmental quality standards mainly take humans and the population as direct risk receptors, and there is no soil environmental quality standard with the groundwater environment as the direct risk receptor. As groundwater is a non-renewable natural resource, it is particularly crucial to assess the soil pollution risks of new pollutants with the groundwater environment as the risk receptor and formulate soil environmental quality standards. Relevant departments are actively promoting the formulation of emission standards and risk control strategies for PPCPs. Antibiotic drugs have been included in the "List of New Pollutants under Key Control", and the "New Pollutant Control Plan" issued in 2022 has successively guided the conduct of PPCPs pollution investigations and ecological risk assessments in multiple basins.

[0021] The "Technical Guidelines for Risk Assessment of Soil Pollution in Construction Land" promulgated by the Ministry of Ecology and Environment in 2019 gives the calculation method for the soil risk control value to protect groundwater as shown in formula (1): (1); Where, is the soil risk control value to protect groundwater, is the maximum concentration limit of the pollutant in groundwater, is the leaching factor for the pollutant in soil to migrate into groundwater.

[0022] The calculation of the leaching factor in formula (1) mainly relies on the soil - water distribution coefficient of the compound, as well as static parameters such as the Darcy rate of groundwater, the thickness of the groundwater mixing zone, the infiltration rate of water in soil, the width of the pollution source area, and the thickness of the underlying polluted soil. The soil risk control value derived from the formula only considers the solid - liquid distribution of pollutants in soil and does not consider the concentration attenuation of soil pollutants during the vertical leaching process in the vadose zone. Researchers used the vadose zone solute transport model HYDRUS - 1D coupled with the groundwater dilution model to construct the soil risk control value based on protecting groundwater at the plot scale. Although the researchers' method considered the transport process of pollutants, due to the lack of spatial differences in the pollutant transport process under different types of soil and hydrogeological environments at the plot scale, the accuracy of the obtained soil risk control value is relatively low, making it difficult to provide a basis for formulating soil environmental quality standards.

[0023] The technical framework of the present invention is as Figure 1 shown, and the main contents include: (1) soil - water distribution coefficient model; (2) vadose zone attenuation model; (3) groundwater dilution model; (4) risk threshold derivation and zoning model.

[0024] (1)Figure 1 The solid-liquid distribution coefficient Kd prediction model in [the above text] is the solid-liquid distribution coefficient prediction model in the previous text. The present invention uses the method of big data analysis of machine learning models to construct a soil-water distribution coefficient prediction model based on the solid-liquid distribution coefficients of different types of pharmaceutical pollutants with a relatively high detection rate in different types of soils and environments. This method can be used for the soil-water distribution coefficients of different types of pharmaceutical pollutants in different types of soils. Compared with the conventional multiple linear model prediction method, the present invention first takes into account the non-linear relationship between the solid-liquid distribution coefficient of the compound and the soil properties and environmental conditions. Secondly, through the application domain analysis, the solid-liquid distribution coefficients of multiple similar compounds can be predicted. Therefore, the uncertainty is lower and the universality is higher.

[0025] (2) Figure 1 The attenuation coefficient AF prediction model in [the above text] is the unsaturated zone attenuation coefficient prediction model in the previous text. In the present invention, when using HYDRUS-1D to simulate the migration process of pollutants in the unsaturated zone, the Python programming language is introduced for assistance, and batch simulation operations are realized through the phydrus library in the Python environment. This method can automatically generate and manage Hydrus input files in batches, efficiently call the Hydrus core program for calculation, and automatically collect and process the output results. Therefore, it effectively improves the calculation efficiency of the HYDRUS-1D single-point model, helps the model of the migration process of multiple points at the regional scale. At the same time, based on the ArcGis platform, the spatial distribution of the attenuation coefficient of the vadose zone at the regional scale is realized.

[0026] (3) Figure 1 The analytical formula of the dilution coefficient DF in [the above text] is the analytical formula of the groundwater dilution coefficient in the previous text. In order to reduce the complexity of the calculation process, the present invention uses the commonly used analytical formula of the groundwater dilution coefficient proposed by the EPA to calculate the dilution coefficient.

[0027] (4)Risk threshold derivation and zoning model. The present invention uses the method of combining the soil pollution risk value and the spatial distribution of the reduction coefficient of the compound in the regional unsaturated zone to conduct the soil pollution risk assessment for protecting groundwater safety. Compared with the conventional groundwater risk assessment, the present invention can provide method support for the formulation of soil environmental quality standards based on groundwater safety and the soil pollution risk assessment of specific regional groundwater safety.

[0028] It should be noted that the present invention does not limit the research area and the types of drugs.

[0029] The following will combine the attached drawings to detail the technical solutions provided by the embodiments of the present invention.

[0030] Figure 2This is a schematic diagram of the process for assessing the pollution risk of pharmaceutical pollutants to groundwater at the regional scale in the present invention, specifically including the following steps: S201: Obtain various soil characteristic parameter data and soil profile data of the research area; the soil characteristic parameter data includes solid-phase physical and chemical properties, adsorption experimental conditions, and molecular descriptors; the soil profile data includes meteorological data and the physical and chemical properties of target PPCPs; the meteorological data includes precipitation, temperature, humidity, and evaporation; the physical and chemical properties of target PPCPs include half-life, Henry's constant, diffusion coefficient, and organic carbon adsorption coefficient.

[0031] Specifically, data collection: Using "soil OR sediment", "sorption", and "PPCPs" as qualifiers, conduct a literature search through the multidisciplinary literature database Web of Sciences. The data included in the present invention needs to meet the following requirements:

[0032] (1) The data needs to come from batch adsorption experiments that comply with the guidelines of the Organization for Economic Co-operation and Development (OECD) or the US Environmental Protection Agency (EPA).

[0033] (2) Soil and sediment samples should be undisturbed soil, clean and intact, without adding or removing any components. (3) The properties of soil and sediment were measured and reported, mainly including pH value, cation exchange capacity, organic carbon content, and clay content, etc.

[0034] (4) Through linear isotherm fitting, the solid-liquid distribution coefficient (Kd) was reported.

[0035] (5) The R2 of the fitted isotherm should be greater than 0.9 to ensure the quality of the adsorption data.

[0036] Data acquisition: With the solid-liquid distribution coefficient as the prediction target, three major categories of a total of 28 kinds of characteristic parameter data were collected. The abbreviations of the characteristic parameter names are shown in Table 1: The first category: solid-phase physical and chemical properties, including soil pH, cation exchange capacity (CEC), organic carbon content (OC), and clay content (CLAY).

[0037] The second category, adsorption experimental conditions, mainly including solution temperature (TEMP) and solid-liquid ratio (RATIO).

[0038] The third category, molecular descriptors, includes Abraham descriptors (E, S, A, B, V), topological descriptor ATSC2dv, octanol-water partition coefficient (logKOW), solubility (logSw), acid dissociation constant (pKa), base dissociation constant (pKb), molecular weight (MW), molecular volume (MV), number of hydrogen bond acceptors (HBA), double bonds (DB), rotatable bonds (RB), number of rings (NR), number of sulfur atoms (NSA), unsaturation index (UI), hydrophilic factor (HF), topological molecular polar surface area (TPSA), molar refractivity (MR), and polarizability (P). Among them, the Abraham descriptors are from the UFZ-LSER database, ATSC2dv is obtained using the Mordred toolkit, and the remaining molecular descriptors are obtained by searching the PubChem, Chemicalize, and CompTox Chemicals Dashboard databases with the name or CAS number of PPCPs. When needed, the organic matter content (OM) is converted to OC using the formula %OM = 1.724×%OC.

[0039] Table 1

[0040] Continued table

[0041] Data preprocessing of soil characteristic parameters is divided into the following six aspects: (1) Unify the units and symbols of each variable, and convert the solid-liquid distribution coefficient to the natural logarithm form (logKd).

[0042] (2) Use the IQR method to find the outliers of logKd and remove the outliers.

[0043] (3) Impute the missing data in the dataset using the RF method.

[0044] (4) Perform multiple collinearity analysis on characteristic variables such as solid-phase properties (pH, CEC, OC, CLAY), experimental conditions (TEMP, RATIO), and molecular descriptors (E, S, A, B, V, ATSC2dv, logKOW, logSW, HBA, DB, RB, NR, NSA, UI, HF, pKa, pKb, MW, MV, TPSA, MR, P) using Kendall correlation analysis.

[0045] (5) Remove the characteristic variables with significant collinear relationships, retain the remaining characteristic variables as input variables, and set logKd as the output variable.

[0046] (6) Divide the dataset into a training set and a test set, where 75% is used as the training set to train the model and adjust the model parameters, and 25% of the dataset is used as the test set to evaluate the model performance.

[0047] Randomly select three sample points for each administrative region / county in the study area through ArcGIS, and then search for soil profile data and meteorological data. After removing the sample points with missing data, 217 sample points are finally retained.

[0048] SoilGrids is a global digital soil mapping system that uses state-of-the-art machine learning algorithms to map the spatial distribution of global soil properties. The SoilGrids prediction model is fitted using more than 230,000 soil profile observations from the WoSIS database and a series of environmental covariates. The covariates are selected from more than 400 environmental layers, which are derived from Earth observation products and other environmental information, including climate, land cover, and topographic morphology. The output of SoilGrids is a global soil property map with a spatial resolution of 250 meters for six standard depth intervals (0−5 cm, 5−15 cm, 15−30 cm, 30−60 cm, 60−100 cm, 100−200 cm).

[0049] In the present invention, the particle size distribution, bulk density, and organic carbon content in six depth intervals of the sample points are extracted through SoilGrids. Since grid data with a depth greater than 200 cm cannot be obtained and the soil properties below 200 cm vary less, it is determined that the data with a depth greater than 200 cm is the same as the data in the 100−200 cm interval.

[0050] The meteorological data mainly includes precipitation, temperature, and relative humidity of 217 sample points from 2012 to 2021. The daily precipitation data uses the TRMM_3B42 product provided by NASA Earhdata with a resolution of 25 km. The monthly average temperature and relative humidity are obtained from a certain environmental data center with resolutions of 1 km and 10 km respectively, and then the monthly potential evapotranspiration is calculated through formula (2). The daily average precipitation and daily average potential evapotranspiration are obtained by averaging the monthly totals. Formula (2) is as follows: (2); where is the monthly potential evapotranspiration in mm, is the monthly average temperature in °C, is the monthly average relative humidity.

[0051] The half-life, Henry's constant, and organic carbon adsorption coefficient of the pollutants were obtained from a certain database. Through the equations provided by existing research, the boiling point, number of atoms, and number of rings of the organic compounds were input to estimate the diffusion coefficient Dw of the pollutants in water and the diffusion coefficient Dg in air.

[0052] S202: Input the data of various soil characteristic parameters into the solid-liquid partitioning coefficient prediction model. Through the solid-liquid partitioning coefficient prediction model, analyze the non-linear relationship between the solid-liquid partitioning coefficient of pharmaceutical pollutants, various soil characteristic parameter data, and environmental conditions to obtain the solid-liquid partitioning coefficient of pharmaceutical pollutants.

[0053] In an exemplary embodiment, the construction process of the solid-liquid partitioning coefficient prediction model specifically includes: obtaining various soil characteristic parameter sample data and dividing the various soil characteristic parameter sample data into a training set and a test set; training the random forest method, extreme gradient boosting method, and support vector regression method respectively through the training set to obtain the hyperparameter combinations of the random forest method, extreme gradient boosting method, and support vector regression method; detecting the prediction accuracy of the random forest method, extreme gradient boosting method, and support vector regression method based on the hyperparameter combinations respectively through the test set; and determining the method with the highest prediction accuracy as the solid-liquid partitioning coefficient prediction model.

[0054] Specifically, the present invention conducts machine learning modeling prediction on the solid-liquid partitioning coefficient of pharmaceutical compounds in soil based on six steps, namely: data search, data acquisition, preprocessing, model construction, model interpretation, and verification with measured data. The method framework diagram of the unsaturated zone attenuation coefficient prediction model provided by the present invention is as Figure 3 shown.

[0055] Model construction: The present invention uses three algorithms with relatively high usage frequencies, namely the random forest method (Random Forest, RF), extreme gradient boosting method (eXtreme Gradient Boosting, XGBoost), and support vector regression method (Support Vector Regression, SVR), to construct a prediction model for the solid-liquid partitioning coefficient of pharmaceutical pollutants. In order to enable the prediction model of the solid-liquid partitioning coefficient of pharmaceutical pollutants to have the optimal prediction effect, its hyperparameters need to be adjusted to the optimal. When tuning the parameters of the prediction model of the solid-liquid partitioning coefficient of pharmaceutical pollutants, first determine the hyperparameter range, then determine the optimal hyperparameters through the grid search method and 9-fold cross-validation method, and finally obtain the hyperparameter combination with the best model prediction effect.

[0056] The machine learning method uses the Python programming language. The RF, XGBoost, and SVR models come from the open-source database Scikit-Learn, and the development environment is Spyder (v5.4.3).

[0057] Model evaluation: After the model training is completed, the test set is used to check the prediction accuracy of the three models, and the coefficient of determination and the mean squared error (MSE) are used to evaluate and obtain the optimal model. Since the test results are related to the random state of data splitting, the random seed is set to 0 (i.e., the data set is randomly divided). The test process is repeated 100 times to repeatedly demonstrate the performance of the two models. R 2 The calculation method of is shown in formula (3): R 2 (3); Among them, is the predicted value of the model, is the true value, is the average value of the true values, R 2 The closer it is to 1, the better the regression fitting curve effect is, N is the total number of sample points, j is the sample point j .

[0058] The calculation method of MSE is shown in formula (4): (4); Among them, is the predicted value of the model, is the true value, n is the total number of simulated sample points, i is the sample point i .

[0059] Model interpretation: In this study, SHapley Additive exPlanations (SHAP) is used for model interpretation to verify the effectiveness of machine learning modeling. The Shapley value uses an additive feature attribution method to explain the contribution of feature variables to the output. The Shapley value is shown in formula (5): (5); Among them, is the Shapley value calculated for the j th row corresponding to the i th feature, is the predicted value of the i th row, is the average value of the predicted values, j is the j th feature, k is the number of feature values. ​

[0060] After calculating the Shapley values of the training set, the importance of feature variables and their impact on the output are proven. The overall contribution of each feature variable is represented by averaging the absolute values of the Shapley values of each feature variable, and the overall contribution of each feature variable is shown in Equation (6): (6); where is the Shapley value corresponding to the i th feature variable of the j th sample, is the average value of the predicted values, is the i th predicted value of the sample, is the contribution of the j th feature variable, n is the number of training samples, k is the number of feature variables.

[0061] Prediction result: The temperature and solid-liquid ratio are set to the most common values in the training set, and an application domain analysis is performed on the data of three PPCPs, namely trimethoprim, nalidixic acid, and niflumic acid, and the surface soil data of 217 sites in the Yangtze River Delta region. The maximum outlier value of the average Euclidean distance calculated for all samples in the solid-liquid distribution coefficient prediction model training set is 1.2303, and the threshold is 1.2303. The dis of the three PPCPs in the prediction samples are shown in Table 2 as the d i , and the results show that the Kds of the three PPCPs in the surface soil of the Yangtze River Delta region can be well predicted by the solid-liquid distribution coefficient prediction model.

[0062] Table 2

[0063] S203: The soil profile data is processed by the unsaturated zone attenuation coefficient prediction model to obtain the unsaturated zone attenuation factors of pharmaceutical pollutants.

[0064] Machine learning modeling is divided into the following steps.

[0065] (1)Data preprocessing Fourteen characteristic variables, such as bulk density, organic carbon content, sand content, clay content, groundwater depth, mean and standard deviation of daily precipitation, mean and standard deviation of daily potential evaporation, organic carbon adsorption coefficient Koc, half-life, Henry's constant, diffusion coefficient Dw in water and diffusion coefficient Dg in air, are used as potential input variables, and AF is transformed into logAF as the output variable. Kendall correlation analysis is used to perform multicollinearity analysis on the 14 potential input variables, and the variables with strong correlation are removed, and the remaining ones are used as input variables. The data set is divided into a training set and a test set, where 75% is used as the training set to train the model and adjust the model parameters, and 25% of the data set is used as the test set to evaluate the model performance.

[0066] (2)Model construction In this study, RF, XGBoost and SVR are used for model construction. After adjusting the hyperparameters to the optimal values, the coefficient of determination ( R 2 )and MSE are used to evaluate the goodness of the model.

[0067] (3)Model interpretation The SHAP method is used to explain the influence of each input variable on the model output.

[0068] Simulate the prediction results.

[0069] In an exemplary embodiment, processing the soil profile data by the unsaturated zone attenuation coefficient prediction model to obtain the unsaturated zone attenuation factor of pharmaceutical pollutants specifically includes: inputting the soil profile data into a one-dimensional hydrological model for simulation to obtain the pollutant concentration in the pore water at a specific depth of the soil; calculating the unsaturated zone attenuation factor of the pharmaceutical pollutants according to the maximum value of the pollutant concentration in the pore water at a specific depth of the soil.

[0070] The calculation method of the unsaturated zone attenuation factor of pharmaceutical pollutants is shown in formula (7): (7); Where, AF is the unsaturated zone attenuation factor of pharmaceutical pollutants, C e is the pollutant concentration in the pore water of the soil surface layer, C d is the maximum value of the pollutant concentration in the pore water at a specific depth of the soil.

[0071] Specifically, the schematic diagram of the research area and sampling points provided by the present invention is as Figure 4As shown, first, soil profile data of the research area are collected, including bulk density, organic carbon content, clay content, sand content, depth, and meteorological data, including precipitation, temperature, humidity, evaporation, and the physical and chemical properties of target PPCPs, including half-life, Henry's constant, diffusion coefficient, organic carbon adsorption coefficient Koc, etc.; subsequently, the above data are input into the HYDRUS-1D model, and the vadose zone attenuation factor AF of pollutants under different environmental conditions is calculated in batches using python; finally, through a machine learning algorithm, with soil profile data, meteorological data, and PPCPs properties as input variables and the unsaturated zone attenuation factor of pharmaceutical pollutants as the output variable, a prediction model for the unsaturated zone attenuation coefficient is constructed.

[0072] HYDRUS-1D is a finite element model used to simulate water, heat, and solute movement in a one-dimensional variably saturated medium. HYDRUS-1D solves the convection-dispersion equation of solute transport and the Richards equation of variably saturated water flow numerically, and considers processes such as root water uptake, infiltration, evaporation, soil water storage, capillary rise, deep drainage, groundwater recharge, and surface runoff. The solute transport equation considers convection and diffusion in the liquid phase, as well as diffusion in the gas phase, and considers linear equilibrium reactions, zero-order generation, and first-order degradation reactions between the liquid and gas phases. HYDRUS-1D also considers different physical and chemical non-equilibrium transport situations, but it is difficult to obtain the corresponding parameters on a regional basis. Therefore, the present invention assumes that the pollutants reach the water-soil distribution equilibrium instantaneously after entering the surface soil, and then enter the aquifer system through the clean vadose zone. The migration process is mainly driven by precipitation convection and hydrodynamic dispersion, without considering physical and chemical non-equilibrium transport.

[0073] In terms of model settings, van Genuchten is selected as the hydraulic model, and Rosetta3 is used to estimate the flow parameters. The upper boundary condition of the flow is set as the atmospheric boundary condition, and the lower boundary condition is set as the free drainage boundary. The upper boundary condition of solute transport is set as the stagnant boundary of volatile solutes, and the lower boundary condition is set as the zero concentration gradient boundary. The upper boundary concentration of the solute is set as 1 mg·L -1 , because the migration process is mainly driven by precipitation, when the daily precipitation is zero, the upper boundary concentration of the solute is set as 0. The unsaturated zone profile is set as 6 m, and the number of profile layers is set as six layers (0−5, 5−15, 15−30, 30−60, 60−100, 100−600 cm). The simulation time is 10 years, and observation points are set at depths of 1 m, 2 m, 3 m, 4 m, 5 m, and 6 m in the profile respectively.

[0074] Use the Phydrus 0.2.0 toolkit in Python to conduct batch HYDRUS-1D simulations for 217 points to obtain the solute concentration at the groundwater level over time. When the pollutant concentration at the groundwater level reaches dynamic equilibrium, take the maximum concentration as the model output. Figure 5 This is the conceptual model diagram for HYDRUS-1D calculation provided by the present invention, as Figure 5 shown. Calculate the unsaturated zone attenuation factor AF through formula (5).

[0075] Conduct HYDRUS-1D simulations on the migration of trimethoprim, nalidixic acid, and niflumic acid in the unsaturated zones of 217 sampling points in the study area to obtain the organic pollutant concentrations at 6 depths at dynamic equilibrium, and calculate the unsaturated zone attenuation factor AF. Subsequently, using the RF, XGBoost, and SVR algorithms, with soil profile data, meteorological data, and organic pollutant properties as input parameters, a machine learning model for predicting the unsaturated zone attenuation factor AF of trimethoprim, nalidixic acid, and niflumic acid was developed.

[0076] The results show that: (1) The AF of deep soil is greater than that of shallow soil, and the AF ranking of the 3 organic pollutants is niflumic acid > nalidixic acid > trimethoprim; (2) The prediction performance of the XGBoost model is better than that of the RF model and the SVR model, and the R2 of the test set reaches 0.9; (3) Through the feature importance analysis based on the SHAP method, it is shown that the profile depth, precipitation, and organic carbon adsorption coefficient are the key factors affecting the unsaturated zone attenuation factor AF (the combined contribution > 76%). In summary, the present invention provides an effective solution for predicting the spatial heterogeneity of AF, and at the same time demonstrates the great potential of combining traditional hydrogeological models with modern machine learning technologies, providing a new method for future research in the environmental field.

[0077] S204: Obtain the groundwater dilution factor of organic pollutants in the study area.

[0078] The calculation method of the groundwater dilution factor of organic pollutants is as shown in formula (8): (8); where, DF is the groundwater dilution factor of organic pollutants, C d is the maximum value of the pollutant concentration in the pore water at a specific depth of the soil, Cgw is the concentration of organic pollutants in the groundwater after dilution.

[0079] C d is the maximum value of the pollutant concentration in the pore water at a specific depth of the soil in the study area, Cgw is the concentration of organic pollutants in the groundwater after dilution in the study area,C d and Cgw Both can be obtained through actual tests or model predictions.

[0080] The calculation method of the dilution factor of organic pollutants in groundwater is shown in formula (9): (9); wherein, K is the hydraulic conductivity of the aquifer, i is the hydraulic gradient, d is the depth of the mixing zone, I is the infiltration rate, L is the length of the source area parallel to the groundwater flow.

[0081] K, i, d, I, L All are parameters corresponding to the study area and can be obtained through actual tests or model predictions.

[0082] S205: Determine the risk value of pharmaceutical pollutants through the solid-liquid partition coefficient of pharmaceutical pollutants, the attenuation factor of pharmaceutical pollutants in the unsaturated zone, and the dilution factor of organic pollutants in groundwater.

[0083] The calculation method of the risk value of pharmaceutical pollutants is shown in formula (10): (10); wherein, is the risk value of pharmaceutical pollutants, is the groundwater concentration limit of pharmaceuticals, K d is the solid-liquid partition coefficient, AF is the attenuation factor in the unsaturated zone, DF is the dilution factor of organic pollutants in groundwater.

[0084] Screening of groundwater concentration limits of three pharmaceutical compounds in groundwater: The LOD of trimethoprim, nalidixic acid, and niflumic acid in groundwater and the PNECwater in the water environment obtained through literature retrieval and formula calculation are shown in Table 3.

[0085] Table 3

[0086] Due to the lack of relevant data for nalidixic acid and niflumic acid, the PNECwater cannot be calculated. Through comparison and in combination with the principle of conservatism, the concentration limits of three PPCPs, namely trimethoprim, nalidixic acid, and niflumic acid, in groundwater are finally determined to be 0.19 ng·L -1 , 1 ng·L -1 and 2 ng·L -1 .

[0087] Derivation results of soil pollution risk values for three pharmaceutical pollutants based on groundwater protection: Since some of the derivation results are too large and lack practical significance, 100 times the highest concentration actually detected of the pollutant in soil or sediment is used as the screening criterion for the derivation result dataset. By excluding extreme values, the derivation dataset is ensured to be closer to the actual situation. After literature research, the maximum concentrations of trimethoprim and nalidixic acid actually detected in soil or sediment are 60 and 455 μg·kg respectively -1 . Since there are few investigations on the concentration of niflumic acid in soil or sediment at present, ibuprofen, which is also a non-steroidal anti-inflammatory drug, is selected as a substitute in this invention, and its maximum detected value is 6.046 mg·kg -1 . After calculation, the final derivation results of the soil safety thresholds based on groundwater protection for the 3 PPCPs are shown in Table 4, and Table 4 is the descriptive statistics of the soil risk values of the 3 pharmaceutical pollutants

[0088] Table 4

[0089] Risk zoning: Risk zoning is carried out according to the attenuation coefficient of the unsaturated zone concentration, and the results are as follows: The attenuation factors log of the 3 pharmaceutical pollutants in the unsaturated zone at different depths AF As Figure 6 shown, when the profile depth is 1 m, the average log AF is 3.08 and the standard deviation is 1.21. When the profile depth is 2 m, the average log AF is 7.30 and the standard deviation is 2.70. When the profile depth is 3 m, the average log AF is 12.10 and the standard deviation is 3.96. When the profile depth is 4 m, the average log AF is 15.50 and the standard deviation is 3.71. When the profile depth is 5 m, the average log AF is 17.99 and the standard deviation is 3.19. When the profile depth is 6 m, the average log AFwas 19.24, and the standard deviation was 2.77. logAF increased significantly with increasing depth (p < 0.05), indicating that the deeper the depth, the lower the concentration of pharmaceutical pollutants. The depth of the unsaturated zone is an important factor affecting the migration and fate of PPCPs in the unsaturated zone. This is consistent with the research results of the researchers. Through field sampling of organochlorine pesticides, it was found that the content of organochlorine pesticides in the soil surface layer was the highest and decreased gradually with the increase of sampling depth. In addition, when the profile depth > 1 m, the error bars of logAF were longer, indicating that the differences in logAF of PPCPs in deep soil were larger, and the migration and fate of PPCPs in deep soil were more easily affected by soil physical and chemical properties and meteorological conditions. The logAF of trimethoprim, nalidixic acid, and niflumic acid was as Figure 7 shown, and the average logAFs were 9.94, 11.09, and 11.56, respectively, and the standard deviations were 6.27, 6.41, and 6.09. This is consistent with the research results of the researchers. In addition to soil conditions and climate conditions, the physical and chemical properties of PPCPs themselves are the key factors affecting the volatilization, degradation, and adsorption of PPCPs in soil. When the half-life is small, the accumulation of PPCPs in soil is high; when the Henry constant is large, the volatility of PPCPs is stronger and the accumulation in soil is low.

[0090] The present invention relates to a method for evaluating the pollution risk of pharmaceutical pollutants to groundwater. The method provided by the present invention is based on big data, coupling a machine learning method and the vadose zone solute transport model HYDRUS-1D, and combines Python software and the ArcGis platform to quantitatively evaluate the groundwater environmental risk and spatial distribution of soil pollution at the regional scale. This method is based on the generalization of the migration process of pollutants in the vadose zone, and quantifies the pollution process of pollutants in the soil-groundwater system, including the solid-liquid distribution coefficient, the concentration attenuation coefficient of the vadose zone, the groundwater dilution factor, etc.; uses the machine learning method to derive the solid-liquid distribution coefficient of pharmaceutical compounds in different types of soil, and uses the HYDRUS-1D model to batch simulate the transport process of pharmaceutical pollutants in the vadose zone under different hydrogeological conditions at the regional scale through the phydrus library in the Python environment, and obtains the spatial distribution of the concentration attenuation coefficient of vadose zone pollutants at the regional scale; at the same time, combines the groundwater environmental quality standard of the target pollutant and hydrogeological parameters such as the regional groundwater depth to evaluate the groundwater environmental risk and spatial zoning of soil pollution at a specific time and region.

[0091] When applying a method for evaluating the pollution risk of pharmaceutical pollutants to groundwater provided by the present invention, it is not necessary to execute according to Figure 2 the order of the steps shown. The specific execution order of each step can be determined according to needs, and the present invention does not limit this.

[0092] The above is a method for assessing the pollution risk of pharmaceutical pollutants to groundwater provided by one or more embodiments of the present invention. Based on the same idea, the present invention also provides a corresponding device for assessing the pollution risk of pharmaceutical pollutants to groundwater, as Figure 8 shown.

[0093] Figure 8 The figure is a schematic diagram of a device for assessing the pollution risk of pharmaceutical pollutants to groundwater at the regional scale provided by the present invention, including: An acquisition module 801, configured to acquire various soil characteristic parameter data and soil profile data of the research area; the soil characteristic parameter data includes solid-phase physical and chemical properties, adsorption experimental conditions, and molecular descriptors; the soil profile data includes meteorological data and the physical and chemical properties of target PPCPs; the meteorological data includes precipitation, temperature, humidity, and evaporation; the physical and chemical properties of target PPCPs include half-life, Henry's constant, diffusion coefficient, and organic carbon adsorption coefficient.

[0094] A first determination module 802, configured to input various soil characteristic parameter data into a solid-liquid partition coefficient prediction model, and analyze the non-linear relationship between the solid-liquid partition coefficient of pharmaceutical pollutants and various soil characteristic parameter data and environmental conditions through the solid-liquid partition coefficient prediction model, so as to obtain the solid-liquid partition coefficient of pharmaceutical pollutants.

[0095] A second determination module 803, configured to process the soil profile data through an unsaturated zone attenuation coefficient prediction model to obtain the unsaturated zone attenuation factor of pharmaceutical pollutants.

[0096] A third determination module 804, configured to obtain the groundwater dilution factor of organic pollutants in the research area.

[0097] A fourth determination module 805, configured to determine the risk value of pharmaceutical pollutants through the solid-liquid partition coefficient of pharmaceutical pollutants, the unsaturated zone attenuation factor of pharmaceutical pollutants, and the groundwater dilution factor of organic pollutants.

[0098] For the specific limitations on the framework of a method for assessing the pollution risk of pharmaceutical pollutants to groundwater at the regional scale, reference can be made to the limitations on the method for assessing the pollution risk of pharmaceutical pollutants to groundwater in the above text, which will not be elaborated here. Each module in the above framework of the method for assessing the pollution risk of pharmaceutical pollutants to groundwater at the regional scale can be implemented in whole or in part through software, hardware, and their combinations. The above modules can be embedded in or independent of the processor in the computer device in the form of hardware, or stored in the memory of the computer device in the form of software, so that the processor can call and execute the operations corresponding to the above modules.

[0099] The present invention also provides a computer-readable storage medium, which stores a computer program, and the computer program can be used to execute the aboveFigure 2 A method for assessing the pollution risk of pharmaceutical pollutants to groundwater at the regional scale is provided.

[0100] The present invention also provides Figure 9 A schematic structural diagram of the computer device as shown, such as Figure 9 shown, at the hardware level, the computer device includes a processor, an internal bus, a network interface, a memory, and a non-volatile memory. Of course, it may also include other hardware required for other services. The processor reads the corresponding computer program from the non-volatile memory into the memory and then runs it to implement the above Figure 2 A method for assessing the pollution risk of pharmaceutical pollutants to groundwater at the regional scale is provided.

[0101] Those of ordinary skill in the art can understand that all or part of the processes in the methods of the above embodiments can be completed by instructing relevant hardware through a computer program. The computer program can be stored in a non-volatile computer-readable storage medium. When the computer program is executed, it can include the processes of the embodiments of the above methods. Among them, any reference to a memory, storage, database, or other medium used in the embodiments provided by the present invention can include at least one of non-volatile and volatile memories. The non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, or optical memory, etc. The volatile memory can include random access memory (RAM) or an external cache memory. By way of illustration and not limitation, RAM can be in various forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM), etc.

[0102] The technical features of the above embodiments can be combined arbitrarily. For the sake of concise description, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, it should be considered to be within the scope recorded by the present invention.

Claims

1. A method for assessing the pollution risk of pharmaceutical pollutants to groundwater at the regional scale, characterized in that, Including: Obtaining various soil characteristic parameter data and soil profile data of the research area; The soil characteristic parameter data includes solid-phase physical and chemical properties, adsorption experimental conditions, and molecular descriptors; the soil profile data includes meteorological data and the physical and chemical properties of target PPCPs; the meteorological data includes precipitation, temperature, humidity, and evaporation; the physical and chemical properties of the target PPCPs include half-life, Henry's constant, diffusion coefficient, and organic carbon adsorption coefficient; Inputting the various soil characteristic parameter data into the solid-liquid partition coefficient prediction model, and analyzing the non-linear relationship between the solid-liquid partition coefficient of pharmaceutical pollutants, various soil characteristic parameter data, and environmental conditions through the solid-liquid partition coefficient prediction model to obtain the solid-liquid partition coefficient of pharmaceutical pollutants; Processing the soil profile data through the unsaturated zone attenuation coefficient prediction model to obtain the unsaturated zone attenuation factor of pharmaceutical pollutants; Obtaining the groundwater dilution factor of organic pollutants in the research area; Determining the risk value of pharmaceutical pollutants through the solid-liquid partition coefficient of pharmaceutical pollutants, the unsaturated zone attenuation factor of pharmaceutical pollutants, and the groundwater dilution factor of organic pollutants.

2. The method according to claim 1, characterized in that, The construction process of the solid-liquid partition coefficient prediction model specifically includes: Obtaining various soil characteristic parameter sample data and dividing the various soil characteristic parameter sample data into a training set and a test set; Training the random forest method, the extreme gradient boosting method, and the support vector regression method respectively through the training set to obtain the hyperparameter combinations of the random forest method, the extreme gradient boosting method, and the support vector regression method; Detecting the prediction accuracy of the random forest method, the extreme gradient boosting method, and the support vector regression method based on the hyperparameter combinations respectively through the test set; Determining the method with the highest prediction accuracy as the solid-liquid partition coefficient prediction model.

3. The method according to claim 1, characterized in that, The processing of the soil profile data through the unsaturated zone attenuation coefficient prediction model to obtain the unsaturated zone attenuation factor of pharmaceutical pollutants specifically includes: Inputting the soil profile data into the one-dimensional hydrological model for simulation to obtain the pollutant concentration in the pore water at a specific depth of the soil; Calculating the unsaturated zone attenuation factor of the pharmaceutical pollutants according to the maximum value of the pollutant concentration in the pore water at a specific depth of the soil.

4. The method according to claim 3, wherein The calculation method of the unsaturated zone attenuation factor of the pharmaceutical pollutants is: ； Among them, AF is the attenuation factor of pharmaceutical pollutants in the unsaturated zone, is the pollutant concentration in the pore water of the soil surface layer, C d is the maximum value of the pollutant concentration in the pore water at a specific depth of the soil.

5. The method according to claim 1, wherein The calculation method of the groundwater dilution factor of organic pollutants is: ； Among them, DF is the dilution factor of organic pollutants in groundwater, Cgw is the concentration of organic pollutants in groundwater after dilution, C d is the maximum value of the pollutant concentration in the pore water at a specific depth of the soil.

6. The method according to claim 1, characterized in that The calculation method of the groundwater dilution factor of organic pollutants is: ； Among them, K is the hydraulic conductivity of the aquifer, i is the hydraulic gradient, d is the depth of the mixing zone, I is the infiltration rate, L is the length of the source area parallel to the groundwater flow.

7. The method according to claim 1, characterized in that The calculation method of the risk value of the pharmaceutical pollutants is: ； Among them, is the risk value of pharmaceutical pollutants, is the groundwater concentration limit of pharmaceuticals, K d is the solid-liquid distribution coefficient of pharmaceutical pollutants, AF is the attenuation factor of pharmaceutical pollutants in the unsaturated zone, DF is the groundwater dilution factor of organic pollutants.

Images