Method for analyzing water pollution in whole chain of river basin based on emission inventory and water quality monitoring

Abstract

This invention discloses a method for analyzing the entire chain of watershed water pollution based on emission inventories and water quality monitoring, belonging to the field of watershed water environment management technology. To address the technical problem of the lack of quantitative link between emission inventories and water quality monitoring, this invention uses emission inventory data, river inflow calculation data, and water quality monitoring data as a foundation. Through multi-level Bayesian optimization and spatial decay modeling, it achieves a quantitative closed-loop analysis method for the entire chain of emission volume → river inflow → monitoring load, providing methodological support for precise watershed pollution control. This invention solves the problem of multi-source coefficient optimization under severely underdetermined conditions at a single cross-section; provides multi-level refined optimization capabilities from global to local levels; reveals the spatial redistribution effect between emission volume and water quality impact; and establishes a complete uncertainty quantification system.

Description

Technical Field

[0001] This invention belongs to the field of watershed water environment management technology, specifically involving a method for analyzing the entire chain of watershed water pollution based on emission inventories and water quality monitoring. Background Technology

[0002] The core objective of watershed water pollution management is to accurately identify pollution sources and establish a quantitative link between discharge volume and water quality response. Currently, discharge inventory compilation technology is relatively mature, and there are studies on gridded discharge inventory accounting in China. However, the existing technology has the following main problems: (1) Disconnection between discharge inventory and monitoring data: Existing discharge inventory studies generally stop at the stage of calculating the amount of discharge into the river, and fail to quantitatively compare and verify the theoretical discharge volume with the measured load of the downstream monitoring section, forming a "broken chain" between the discharge inventory and water quality monitoring. (2) Lack of data-driven verification for the discharge coefficient: The discharge coefficient usually adopts the empirical value within the range recommended by the technical guidelines, but in actual application, due to problems such as improper coefficient values ​​and data entry errors, there is often an order of magnitude deviation between the theoretical discharge volume and the actual monitoring value. (3) Failure to consider the spatial attenuation effect during river transport: The attenuation rate of different pollutants in the river is significantly different. Ignoring the spatial distance attenuation will lead to a systematic overestimation of the impact of distant pollution sources on the water quality of the monitoring section. (4) Coefficient optimization is a severely underdetermined problem: watersheds usually have only a limited number of monitoring sections (few constraint equations), while the number of parameters to be optimized (multiple inflow coefficients of multiple pollution sources) is far greater than the number of constraint equations, and traditional optimization methods are difficult to provide robust solutions. (5) The internal parameters of complex pollution sources such as urban non-point sources cannot be optimized in a refined manner: global optimization only provides an overall correction factor and cannot distinguish the contribution differences of different functional areas and the sources of parameter deviations. Summary of the Invention

[0003] To address the technical challenge of the lack of a quantitative link between emission inventories and water quality monitoring, this invention aims to provide a method for analyzing the entire chain of watershed water pollution based on emission inventories and water quality monitoring. Using emission inventory data, river inflow calculation data, and water quality monitoring data as a foundation, this method employs multi-level Bayesian optimization and spatial decay modeling to achieve a quantitative closed-loop analysis of the entire chain from "emissions to river inflows to monitoring load," providing methodological support for precise watershed pollution control.

[0004] This invention provides a method for analyzing the entire chain of watershed water pollution based on emission inventories and water quality monitoring, comprising the following steps:

[0005] Step 1: Estimation of discharge and inflow into rivers

[0006] (a) Emissions estimation

[0007] For point sources such as industrial enterprises and urban domestic sewage treatment facilities, pollutant emissions are calculated based on online monitoring data of the pollution source, discharge permit data, or environmental statistics.

[0008] For agricultural and rural sources (including free-range livestock and poultry farming, aquaculture, agricultural planting, and rural life) and domestic sources (urban domestic scattered emissions) in non-point sources, as well as large-scale livestock and poultry farming in point sources, the production and emission coefficient method is adopted, which is calculated based on pollution source activity level data, sewage collection and treatment data, and corresponding production and emission coefficients.

[0009] For urban non-point source pollution (urban runoff pollution), a simplified model method is used to calculate the total amount of pollutants discharged during rainfall in different land types of catchment areas based on the catchment area, annual rainfall, runoff coefficient and average pollutant concentration of different functional areas of the city.

[0010] (b) Estimation of pollutant discharge into the river: Based on the pollutant discharge amount calculated by source, the discharge correction coefficient for each type of pollution source is determined according to factors such as the distance between the pollution source discharge outlet and the receiving water body, the degree of pipelineization of the transmission route, regional topography and climate conditions, and the discharge amount is multiplied by the discharge correction coefficient to obtain the pollutant discharge into the river.

[0011] Step 2: Intelligent preprocessing and missing data completion for monitoring data

[0012] (a) Water quality index filling: Linear interpolation was used to fill hourly missing data for chemical oxygen demand (COD), ammonia nitrogen (NH3-N), total nitrogen (TN), and total phosphorus (TP);

[0013] (b) Three-level intelligent data filling for traffic flow:

[0014] Level 1: Fill in the data by matching the median traffic for the same month, day, and time from historical data over many years;

[0015] Level 2: If there is no matching data in Level 1, the median traffic of the same day and time in adjacent months will be used.

[0016] Level 3: If there is no match in the first two levels, linear interpolation is used;

[0017] (c) Outlier detection: A threshold method based on physical meaning (e.g., 0~50 m³ / s) is used instead of the statistical IQR method to preserve the true flood peak events;

[0018] (d) Missing data sensitivity analysis: Three scenario analysis frameworks were established: conservative scenario (missing data = 0), medium scenario (missing data = average hourly load), and high scenario (missing data = 1.5 times the average hourly load). The impact of missing data on load estimation was quantified, and the medium scenario was selected as the subsequent optimization target.

[0019] Step 3: Monitoring load calculation and establishment of the whole-chain coefficient system

[0020] (a) Based on hourly monitoring data, the annual pollutant load is calculated using the integral method:

[0021]

[0022] Where L is the annual pollutant load, kg; C i Q represents the pollutant concentration monitored in the i-th hour, in mg / L; i Let m be the instantaneous flow rate monitored in the i-th hour. 3 / s; i ranges from 1 to n, where n is the number of valid monitoring hours throughout the year, and 3.6 is the unit conversion factor (from 3600s / h × 10). -3 The value is derived from kg / mg and is used to convert mg·m³ / (L·s) to kg.

[0023] (b) Establish a three-level coefficient system of “discharge → inflow → monitoring load”, and calculate the comprehensive inflow coefficient (inflow / discharge), river transport coefficient (monitoring load / inflow) and comprehensive transport coefficient (monitoring load / discharge) respectively.

[0024] Step 4: Bayesian Global Coefficient Optimization

[0025] For the severely underdetermined problem (m pollution sources, only one constraint equation for each pollutant, and m+1 parameters to be estimated), a Bayesian maximum a posteriori (MAP) estimation framework is adopted:

[0026] (a) Problem formalization:

[0027]

[0028] Among them, M p To supplement the load of the monitoring section, Let p be the emission amount of pollutant from the j-th source. Let p be the initial inflow coefficient of pollutant p from the j-th source. U is the correction factor for the inflow coefficient of pollutant p from the j-th source to be optimized. p Contribute to unknown sources that were not included in the list; j takes values ​​from 1 to m.

[0029] (b) Differentiated prior distribution:

[0030] The correction factor uses a truncated normal prior, in the following form:

[0031]

[0032] Prior mean μ of each source j and variance σ j2 The standard deviation is set according to the data reliability (a small standard deviation is used for sources with high data reliability, and a large standard deviation is used for sources with high uncertainty); a and b are the lower and upper limits of the truncated normal distribution, which are used to constrain the reasonable physical range of the correction factor.

[0033] Unknown source U p Using Gamma distribution prior The shape parameter ν and the rate parameter ω are set to encourage smaller contributions from unknown sources.

[0034] (c) Likelihood function: Assuming the observation error follows a normal distribution, the likelihood function is:

[0035]

[0036] In the formula: For parameter-based The model predicts the load; Let M be the standard deviation of the observation error. p The assumption is 10%. This reasonably covers the uncertainties inherent in the monitoring data itself, including instrument errors, sampling representativeness, and flow measurement errors.

[0037] (d) Solution algorithm: The differential evolution global optimization algorithm is used to solve the MAP estimation, and the 95% confidence interval of the parameters is calculated by the Laplace approximation (the inverse of the Hessian matrix);

[0038] (e) Anomaly source detection: Calculate the correction factor for each source.

[0039]

[0040] In the formula: μ is the optimal estimate of the j-th source correction factor. j and σ j Let z be the mean and standard deviation of its prior distribution. j When the threshold is exceeded (e.g., 2), the source is marked as an abnormal source, providing a priority ranking for data verification.

[0041] Step 5: Secondary optimization of functional zones for complex pollution sources

[0042] For pollution sources with complex internal structures (such as urban non-point sources), the global correction factor is decomposed to the functional zone level:

[0043] (a) Model establishment: The formula for calculating urban non-point source pollution inflow into rivers is:

[0044]

[0045] Among them, W a The effective weight of the a-th functional area, , where p is the event-means concentration (EMC) of pollutant p in functional zone a, in mg / L; a ranges from 1 to h, where h is the total number of functional zones; R is the stormwater pipe network correction factor, and D is the distance correction factor. The initial values ​​of R and D can be determined based on the "Technical Guidelines for Compiling Water Pollution Source Discharge Inventory" or local experience, and are used as parameters to be estimated in the optimization.

[0046] (b) Parameter optimization: There are a total of 4h+2 parameters (cj of h functional zones × 4 pollutants, plus R and D), 4 constraint equations (target inflow of each pollutant into the river), which are solved using the same Bayesian MAP framework as in step four, with the prior distribution center set as the original coefficient value;

[0047] (c) Contribution structure redistribution: Analyze the changes in the contribution of each functional zone to the river inflow before and after optimization, and identify the main sources of parameter deviation.

[0048] Step Six: Monte Carlo Uncertainty Analysis

[0049] (a) Set the emissions from each pollution source to follow a normal distribution (coefficient of variation CV), and the inflow coefficient to follow a uniform distribution within a multiple range, and conduct N (e.g., 10,000) random sampling simulations;

[0050] (b) Obtain the probability distribution characteristics (mean, standard deviation, quantiles) of the total amount of each pollutant entering the river, and compare it with the monitoring load to assess the degree of closure;

[0051] (c) Perform elasticity coefficient analysis: Increase the emission amount by a fixed percentage (e.g., 10%) at each source and calculate the proportion of the impact on the total amount entering the river. Identify the most sensitive sources of pollution.

[0052] Step 7: Establishment of Spatial Distance-Attenuation Model and Analysis of Effective Contribution

[0053] (a) Spatial decay model:

[0054]

[0055] In the formula: To supplement the load of the monitoring section, Let p be the emission amount of pollutant from the j-th source. Let k be the initial inflow coefficient of pollutant p from the j-th source. p The river channel decay rate of pollutant p, in km -1 ;d j γ is the river channel distance from the j-th source to the outflow monitoring section, in km; pis the global correction factor for pollutant p; for area sources, the representative distance is taken by aggregating according to the control unit, and for point sources, the distance to the monitoring section is calculated for each individual source; the value of j is 1~m;

[0056] (b) Joint optimization: The differential evolution algorithm is used to simultaneously optimize the decay rate k. p and global correction factor γ p This allows the model predictions to approximate the monitoring load after completion; each pollutant is optimized independently.

[0057] (c) Effective contribution rearrangement: After considering the spatial distance attenuation, the effective contribution rate of each source to the monitoring section is recalculated to reveal the spatial redistribution effect between "emission contribution" and "water quality impact".

[0058] The beneficial effects of this invention are:

[0059] (1) For the first time, a complete quantitative closure of the emission inventory and water quality monitoring chain has been achieved. A complete quantitative chain of "emission amount → river inflow amount → monitoring load" has been established, filling the technical gap that existing emission inventory research stops at the calculation of river inflow amount and cannot be connected with water quality monitoring.

[0060] (2) The problem of multi-source coefficient optimization under severe underdetermined conditions of a single section is solved. By integrating differentiated prior knowledge and observation constraints through a Bayesian framework, a physically reasonable and quantifiable correction factor can still be obtained under the extreme underdetermined condition of 1 constraint equation / m+1 unknowns.

[0061] (3) It provides multi-level refined optimization capabilities from global to local. The second-level Bayesian optimization decomposes the global correction factor to the functional area level (4h+2 parameters), revealing the differences in the internal structure of the river inflow coefficient and the specific sources of parameter deviation, realizing refined diagnosis from "which type of source needs to be corrected" to "which parameter inside the source needs to be corrected".

[0062] (4) It reveals the spatial redistribution effect between emissions and water quality impact. The spatial distance-attenuation model quantitatively confirms the differentiated attenuation characteristics of different pollutants, reveals the significant separation between "major emitters" and "major water quality impactors", and provides a scientific basis for differentiated governance strategies.

[0063] (5) A complete uncertainty quantification system was established. Monte Carlo simulation (10,000 times) provides the probability distribution characteristics of the amount of water entering the river, elasticity coefficient analysis identifies the most sensitive pollution sources, missing sensitivity analysis quantifies the impact of data incompleteness, z-score anomaly detection marks suspicious data sources, and multiple independent methods are used to cross-validate the reliability of the analysis results.

[0064] (6) The three-level data quality verification method effectively identifies hidden errors in the river inflow coefficient data. Through compliance checks, cross-table verification, and outlier investigation, problems such as classification errors, column misalignment, and duplicate entries that are difficult to identify by conventional checks can be found. Attached Figure Description

[0065] Figure 1 This is the overall flowchart of the full-chain closure analysis method of the present invention;

[0066] Figure 2 A schematic diagram illustrating the hierarchical relationship between Bayesian global coefficient optimization and secondary functional area optimization;

[0067] Figure 3 This is a schematic diagram illustrating the working principle of the spatial distance-attenuation model.

[0068] Figure 4 The Bayesian correction factor heatmap for Example 1 (9 types of pollution sources × 4 types of pollutants);

[0069] Figure 5 This is a comparison diagram of the spatial effective contribution rearrangement in Example 1. Detailed Implementation

[0070] The present invention will be further illustrated by the following embodiments, but is not limited to the following embodiments.

[0071] Example 1:

[0072] Taking the Nanchuan River Basin in Zhongyang County, Shanxi Province as an example, the specific implementation process of the method of the present invention is explained.

[0073] Basin Overview: The Nanchuan River Basin is located in the hilly and gully region of the Loess Plateau in the middle reaches of the Yellow River, with a basin area of ​​1,438.62 km². 2 It involves 6 townships, and the watershed can be divided into 4 control units.

[0074] Data source:

[0075] Emission inventory data: Basic data for industrial sources comes from the key monitoring unit platform and Shanxi Provincial Environmental Statistics Data; basic data for urban domestic sources (centralized sewage treatment facilities and decentralized urban domestic emissions) comes from Shanxi Provincial Environmental Statistics Data, online monitoring data of urban sewage treatment plants, and the Lüliang City Environmental Statistics Yearbook; basic data for large-scale livestock and poultry farming and livestock and poultry farmers (free-range livestock and poultry) comes from on-site survey data of local ecological and environmental departments; basic data for agricultural planting comes from the third land use type survey data of Zhongyang County and the Lüliang City Environmental Statistics Yearbook; data for aquaculture comes from the aquaculture activity level survey data provided by the agricultural department; data for rural domestic sources comes from rural environmental remediation data and the seventh national population census data; and data for urban non-point source pollution comes from the third land use type survey data, the Lüliang City Water Resources Bulletin, and runoff concentration parameters from relevant literature.

[0076] River discharge data: River discharge data for 9 types of pollution sources (calculated based on discharge amount × discharge coefficient × correction coefficient), including four pollutants: COD, NH3-N, TN, and TP;

[0077] Water quality monitoring data: Hourly monitoring data from automatic monitoring stations at the outflow section (2022, COD, NH3-N, TN, TP concentrations and instantaneous flow rates, a total of 5,928 valid records in 2022).

[0078] The overall flowchart of the full-chain closure analysis method of this invention is shown below. Figure 1 As shown, the specific steps include:

[0079] Step 1: Estimate the discharge and inflow volumes.

[0080] (a) Emissions estimation

[0081] For point sources, industrial sources and urban domestic sewage treatment facilities, pollutant emissions are calculated based on online monitoring data of the pollution sources, discharge permit data, or environmental statistics.

[0082] For agricultural and rural sources and domestic sources (urban domestic decentralized emissions) in non-point sources, and large-scale livestock and poultry farming in point sources, the production and emission coefficient method is adopted, which is calculated based on pollution source activity level data, sewage collection and treatment data and corresponding production and emission coefficients.

[0083] For urban non-point source pollution, a simplified model method is used to calculate the total amount of pollutants discharged during rainfall in different land types' catchment areas based on the catchment area, annual rainfall, runoff coefficient, and average pollutant concentration of different functional zones in the city.

[0084] (b) Estimation of pollutant discharge into the river: Based on the source-specific calculation of pollutant discharge, the discharge correction coefficient for each type of pollution source is determined according to the distance between the pollution source discharge outlet and the receiving water body, the degree of pipelineization of the transmission route, regional topography and climate conditions. The discharge amount is multiplied by the discharge correction coefficient to obtain the pollutant discharge into the river.

[0085] Taking the Nanchuan River Basin as an example, based on multi-source data such as key monitoring unit platforms, environmental statistical yearbooks, and on-site surveys, the emissions and discharges into the river from nine types of pollution sources (industrial sources, urban domestic sewage treatment facilities, large-scale livestock and poultry farming, agricultural planting, rural domestic activities, aquaculture, free-range livestock and poultry, decentralized urban domestic emissions, and urban non-point source pollution) were calculated. Figure 2 (Framework of the Full-Chain Quantitative Analysis Technical Route) As shown in the "Data Input" section on the left, this step integrates emission inventories, river discharge calculations, and GIS spatial data, forming the data foundation for subsequent analysis. The calculation results are as follows: COD emission load 1036.52 tons, river discharge load 463.89 tons; NH3-N emission 19.95 tons, river discharge 6.84 tons; TN emission 131.30 tons, river discharge 69.15 tons; TP emission 13.96 tons, river discharge 5.32 tons. These results provide a starting point for establishing a full-chain quantitative relationship of "emission → river discharge → monitoring".

[0086] Step 2: Preprocess the monitoring data.

[0087] (a) Water quality index filling: Linear interpolation was used to fill hourly missing data for chemical oxygen demand (COD), ammonia nitrogen (NH3-N), total nitrogen (TN), and total phosphorus (TP);

[0088] (b) Flow data are filled using three levels of intelligent filling (historical median for the same period → median for adjacent months → linear interpolation), while water quality data (COD, NH3-N, TN, TP) are filled using linear interpolation.

[0089] (c) Outlier detection and missing value sensitivity analysis: Outliers were detected using a threshold method based on physical meaning (0~50 m). 3 / s) detection to preserve the true flood peak. A four-scenario analysis framework was established: conservative, moderate, high, and flood season weighted, to quantify the impact of missing data on load estimation. Conservative scenario: The load for the missing period is zero, and the annual load is the cumulative measured value of the period with data. Moderate scenario: It is assumed that the average hourly load for the missing period is the same as that for the period with data, and each month is amplified independently. Flood season weighted scenario: The flood season and non-flood season are distinguished, and the missing periods are filled in with seasonal weights. First, the average hourly load of the non-flood season (January to June) is calculated as the baseline. Second, for each non-flood season month (January to June and December), the missing periods are filled with the average hourly load of that month. Finally, for each flood season month (July to November), due to the extremely small amount of observational data and the possibility that it is biased towards the low load period when equipment is working normally (high equipment failure rate during rainstorms), the missing periods are filled with the average hourly load of the non-flood season. High scenario: The load for the missing period is 1.5 times the average hourly load of the period with data.

[0090] Step 3: Calculate monitoring load and establish a full-chain coefficient system.

[0091] (a) Based on hourly monitoring data, the annual pollutant load is calculated using the integral method:

[0092]

[0093] Where L is the annual pollutant load, kg; C i Q represents the pollutant concentration monitored in the i-th hour, in mg / L; i Let m be the instantaneous flow rate monitored in the i-th hour. 3 / s; i ranges from 1 to n, where n is the number of valid monitoring hours throughout the year, and 3.6 is the unit conversion factor (from 3600s / h × 10). -3 The value is derived from kg / mg and is used to convert mg·m³ / (L·s) to kg.

[0094] (b) Establish a three-level coefficient system of “discharge → inflow → monitoring load”, and calculate the comprehensive inflow coefficient (inflow / discharge), river transport coefficient (monitoring load / inflow) and comprehensive transport coefficient (monitoring load / discharge) respectively.

[0095] Based on the hourly data supplemented in step two, the annual pollutant load was calculated using the integral method, and a three-level coefficient system of "emissions → river discharge → monitoring load" was established. Finally, the supplemented monitoring load under a moderate scenario was selected as the subsequent optimization target: COD 164.26 tons, NH3-N 5.73 tons, TN 72.49 tons, TP 0.84 tons. Figure 2 As shown, the outputs of this step, "monitoring load target" and "full-chain coefficient system", are directly used as inputs and constraints for subsequent Bayesian optimization.

[0096] The combined inflow coefficients are COD 0.448, NH3-N 0.342, TN 0.526, and TP 0.381; the channel transport coefficients are COD 0.354, NH3-N 0.839, TN 1.049, and TP 0.158.

[0097] Step 4: Bayesian global optimization.

[0098] For the severely underdetermined problem (m pollution sources, only one constraint equation for each pollutant, and m+1 parameters to be estimated), a Bayesian maximum a posteriori (MAP) estimation framework is adopted:

[0099] (a) Problem formalization:

[0100]

[0101] Among them, M p To supplement the load of the monitoring section, Let p be the emission amount of pollutant from the j-th source. Let p be the initial inflow coefficient of pollutant p from the j-th source. U is the correction factor for the inflow coefficient of pollutant p from the j-th source to be optimized. p Contribute to unknown sources that were not included in the list; j ranges from 1 to m;

[0102] (b) Differentiated prior distribution:

[0103] The correction factor uses a truncated normal prior, in the following form:

[0104]

[0105] Prior mean μ of each source j and variance σ j 2 The standard deviation is set according to the data reliability (a small standard deviation is used for sources with high data reliability, and a large standard deviation is used for sources with high uncertainty); a and b are the lower and upper limits of the truncated normal distribution, which are used to constrain the reasonable physical range of the correction factor.

[0106] Unknown source U p Using Gamma distribution prior The shape parameter ν and the rate parameter ω are set to encourage smaller contributions from unknown sources.

[0107] (c) Likelihood function: Assuming the observation error follows a normal distribution, the likelihood function is:

[0108]

[0109] In the formula: For parameter-based The model predicts the load; Let M be the standard deviation of the observation error. p The assumption is 10%. This reasonably covers the uncertainties inherent in the monitoring data itself, including instrument errors, sampling representativeness, and flow measurement errors.

[0110] (d) Solution algorithm: The differential evolution global optimization algorithm is used to solve the MAP estimation, and the 95% confidence interval of the parameters is calculated by the Laplace approximation (the inverse of the Hessian matrix);

[0111] (e) Anomaly source detection: Calculate the correction factor for each source.

[0112]

[0113] In the formula: μ is the optimal estimate of the j-th source correction factor. j and σ j Let z be the mean and standard deviation of its prior distribution. j When the threshold is exceeded (e.g., 2), the source is marked as an abnormal source, providing a priority ranking for data verification.

[0114] This step aims to address the "underdetermined" problem in optimizing multi-source parameters under single-section constraints. For example... Figure 2 As shown, with monitoring load as the target and the emissions from each source and the initial inflow coefficient as inputs, under differentiated truncated normal prior distributions (e.g., prior mean 0.8, standard deviation 0.3 for large-scale livestock and poultry farming; prior mean 1.0, standard deviation 0.3 for rural life), the differential evolution algorithm (seed 42, maximum iterations 1000) is used to solve for the maximum a posteriori estimate, obtaining the correction factor for each source. The optimization results of the correction factor for each pollutant from each pollution source are shown in [the table / example]. Figure 4 A total of 36 correction factors were used for 4 pollutants from 9 sources. *Marked areas in the figure indicate significantly adjusted data. The correction factors for COD and TP from large-scale livestock and poultry farming were significantly reduced to approximately 0.1. The z-score calculation results show that the z-scores for COD and TP from large-scale livestock and poultry farming were 2.33, exceeding the threshold of 2, and were therefore marked as anomalous sources. After optimization, the deviation between the model-predicted load and the monitored load was controlled within 12% (COD +8.8%, NH3-N +3.4%, TN -0.7%, TP +11.4%), validating the effectiveness of the Bayesian framework in optimizing coefficients under underdetermined conditions.

[0115] Step 5: Perform secondary optimization of urban surface sources by functional zone.

[0116] For urban non-point source pollution sources with complex internal structures, the global correction factor is decomposed to the functional zone level:

[0117] (a) Model establishment: The formula for calculating urban non-point source pollution inflow into rivers is:

[0118]

[0119] Among them, W a The effective weight of the a-th functional area, , where p is the event-means concentration (EMC) of pollutant p in functional zone a, in mg / L; h is the total number of functional zones; R is the stormwater network correction factor, and D is the distance correction factor. The initial values ​​of R and D can be determined based on the "Technical Guidelines for the Compilation of Water Pollution Source Discharge Inventory" or local experience, and are used as parameters to be estimated in the optimization.

[0120] (b) Parameter optimization: There are a total of 4h+2 parameters (cj of h functional zones × 4 pollutants, plus R and D), 4 constraint equations (target inflow of each pollutant into the river), which are solved using the same Bayesian MAP framework as in step four, with the prior distribution center set as the original coefficient value;

[0121] (c) Contribution structure redistribution: Analyze the changes in the contribution of each functional zone to the river inflow before and after optimization, and identify the main sources of parameter deviation.

[0122] like Figure 2 As shown, this step extracts the sources to be refined (such as urban non-point source pollution) from the global optimization results in step four, and performs functional zone decomposition and parameter set reconstruction. A two-level optimization model (18 parameters in total) is established, containing 4 (functional zones: commercial zone, industrial zone, residential zone, and transportation zone) × 4 (pollutant) runoff concentration parameters, 1 stormwater network coefficient, and 1 distance coefficient, and is solved again under the Bayesian MAP framework. After optimization, the effective weights of each functional zone (commercial zone, industrial zone, residential zone, and transportation zone) in terms of river inflow are 5.1%, 70.6%, 10.7%, and 13.7%, respectively. After optimization, the COD Cj in the industrial zone decreased from 139 to 34.4 mg / L (-75.3%), TP Cj decreased from 0.80 to 0.26 mg / L (-67.7%), the stormwater network coefficient decreased from 0.80 to 0.74, and the distance coefficient decreased from 0.90 to 0.85.

[0123] Step Six: Monte Carlo Uncertainty Analysis

[0124] (a) Set the emissions from each pollution source to follow a normal distribution (coefficient of variation CV), and the inflow coefficient to follow a uniform distribution within a multiple range, and conduct N (e.g., 10,000) random sampling simulations;

[0125] (b) Obtain the probability distribution characteristics (mean, standard deviation, quantiles) of the total amount of each pollutant entering the river, and compare it with the monitoring load to assess the degree of closure;

[0126] (c) Perform elasticity coefficient analysis: Increase the emission amount by a fixed percentage (e.g., 10%) at each source and calculate the proportion of the impact on the total amount entering the river. Identify the most sensitive sources of pollution.

[0127] like Figure 2 As shown in "Step 4: Monte Carlo Uncertainty Analysis" of the flowchart, this step is independent of Bayesian optimization and evaluates the overall uncertainty of the river inflow accounting system through forward simulation. It is assumed that the emissions from each source follow a normal distribution (CV=20%), and the river inflow coefficient is uniformly distributed within ±50%, and 10,000 random sampling simulations are performed. Simulation results show that the total amount of COD and TP entering the river is greater than the monitored values ​​in 100% of the simulations, indicating a systematic overestimation; NH3-N has the best closure (82.4% greater than the monitored values); and TN has the most symmetrical distribution (39.1% greater than the monitored values). Elasticity coefficient analysis (see...) Figure 2 The "source-by-source sensitivity ranking" results show that large-scale livestock and poultry farming is the largest source of uncertainty in the total amount of COD (39.1%) and TP (52.8%) entering the river, which is consistent with the results of the abnormal sources identified in step four by Bayesian optimization.

[0128] Step 7: Establishment of Spatial Distance-Attenuation Model and Analysis of Effective Contribution

[0129] (a) Spatial decay model:

[0130]

[0131] Where: M p To supplement the load of the monitoring section, Let p be the emission amount of pollutant from the j-th source. Let k be the initial inflow coefficient of pollutant p from the j-th source. p The river channel decay rate of pollutant p, in km -1 ;d j γ is the river channel distance from the j-th source to the outflow monitoring section, in km; p The global correction factor for pollutant p is used; for area sources, the representative distance is obtained by aggregating them according to the control unit, and for point sources, the distance to the monitoring section is calculated for each individual source.

[0132] (b) Joint optimization: The differential evolution algorithm is used to simultaneously optimize the decay rate k. p and global correction factor γ p This allows the model predictions to approximate the monitoring load after completion; each pollutant is optimized independently.

[0133] (c) Effective contribution rearrangement: After considering the spatial distance attenuation, the effective contribution rate of each source to the monitoring section is recalculated to reveal the spatial redistribution effect between "emission contribution" and "water quality impact".

[0134] like Figure 3 As shown, the model quantifies the contribution of pollution sources (point and area source control units) to the monitoring section after attenuation over river distance. Area sources are aggregated into four control units (representing distances of 8, 12, 22, and 30 km), while the river distance to the monitoring section for each point source is calculated individually (range 2.3–99.6 km). Joint optimization yields the following attenuation parameters: TP half-life distance 3.1 km, COD 4.8 km, NH3-N 9.3 km, and TN 11.9 km. The core output of the model is the redistribution of the "effective contribution" of pollution sources, as shown in the figure. Figure 5 (Comparison chart of effective contribution percentage by source) As shown in the chart, the "emission percentage" (black bar ■) and "effective contribution percentage" (white bar □) are compared using a double-column stacked bar chart. It can be clearly seen that for COD, large-scale livestock and poultry farming, which accounts for 56.6% of emissions, has an effective contribution of only 10.3% due to its greater spatial distance; while centralized facilities, which account for only 6.8% of emissions, have an effective contribution as high as 26.8% due to their closer proximity. For example... Figure 5 The illustration visually reveals the significant separation between "major emitters" and "major contributors to water quality."

[0135] The effectiveness of this method was verified in the following ways:

[0136] (1) The deviation between the Bayesian optimized prediction value and the monitored load is controlled within 12%;

[0137] (2) The spatial decay model fitting deviation is 0.00% (perfect fit);

[0138] (3) The probability distribution of the Monte Carlo simulation and the location relationship of the monitored values ​​are consistent with the physical expectations;

[0139] (4) Weak prior sensitivity analysis shows that the anomaly detection results are determined by data rather than prior assumptions.

Claims

1. A method for analyzing the entire chain of watershed water pollution based on emission inventories and water quality monitoring, characterized in that... Includes the following steps: Step 1: Calculate the amount of water pollutants discharged and discharged into rivers in the study area; Step 2 involves intelligent preprocessing of the monitoring data. The traffic data is filled using a three-level intelligent filling strategy based on the median of the same period over many years, and a multi-scenario sensitivity analysis framework is established to determine the optimization target. Step 3: Calculate the annual pollutant load based on time-period monitoring data, and establish a full-chain coefficient system between emissions, river discharge, and monitoring load; Step 4: Using the Bayesian maximum a posteriori estimation framework, the differential prior distribution of each pollution source and the monitoring load constraint are integrated to globally optimize the correction factor of the river inflow coefficient of each source, and suspicious data sources are marked by z-score anomaly detection. Step 5: Perform secondary Bayesian optimization of pollution sources with complex internal structures by functional zones, decomposing the global correction factor into specific parameters at the functional zone level; Step 6: Quantify the uncertainty of the estimated inflow volume through Monte Carlo simulation and elasticity coefficient analysis, and identify the most sensitive pollution sources; Step 7: Construct a spatial distance-attenuation model, jointly optimize the attenuation rate and global correction factor, and calculate the effective contribution rate of each source to the monitoring section.

2. The method for analyzing the entire chain of watershed water pollution based on emission inventories and water quality monitoring according to claim 1, characterized in that, In the first step of emission estimation, point sources and non-point sources are included. Point sources include industrial enterprises, urban sewage treatment facilities, and large-scale livestock and poultry farming. Agricultural and rural sources in non-point sources include free-range livestock and poultry farming, aquaculture, agricultural planting, and rural life. Domestic sources in non-point sources are the decentralized emissions from urban life. Urban non-point sources are urban runoff pollution.

3. The method for analyzing the entire chain of watershed water pollution based on emission inventories and water quality monitoring according to claim 2, characterized in that, In step one, the specific method for calculating the discharge and river discharge is as follows: (a) Emissions estimation For point sources such as industrial enterprises and urban domestic sewage treatment facilities, pollutant emissions are calculated based on online monitoring data of the pollution source, discharge permit data, or environmental statistics. For agricultural and rural sources and domestic sources in non-point sources, and large-scale livestock and poultry farming in point sources, the production and discharge coefficient method is adopted, which is calculated based on pollution source activity level data, sewage collection and treatment data and corresponding production and discharge coefficients. For urban non-point source pollution, a simplified model method is used to calculate the total amount of pollutants discharged during rainfall in different land types' catchment areas based on the catchment area, annual rainfall, runoff coefficient, and average pollutant concentration of different functional zones in the city. (b) Estimation of pollutant discharge into the river: Based on the pollutant discharge amount calculated by source, the discharge correction coefficient for each type of pollution source is determined according to the distance between the discharge outlet of the pollution source and the receiving water body, the degree of pipelineization of the transmission route, regional topography and climate conditions. The discharge amount is multiplied by the discharge correction coefficient to obtain the pollutant discharge into the river.

4. The method for analyzing the entire chain of watershed water pollution based on emission inventories and water quality monitoring according to claim 1, characterized in that, In step two, the specific methods for intelligent preprocessing and missing data completion of monitoring data are as follows: (a) Water quality index filling: Linear interpolation was used to fill hourly missing data for chemical oxygen demand, ammonia nitrogen, total nitrogen, and total phosphorus concentrations; (b) Three-level intelligent data filling for traffic flow: Level 1: Fill in the data by matching the median traffic for the same month, day, and time from historical data over many years; Level 2: If there is no matching data in Level 1, the median traffic of the same day and time in adjacent months will be used. Level 3: If there is no match in the first two levels, linear interpolation is used; (c) Outlier detection: A threshold method based on physical meaning is used instead of the statistical IQR method to preserve real flood peak events; (d) Missing data sensitivity analysis: Three scenario analysis frameworks are established: conservative scenario, medium scenario and high scenario. The impact of missing data on load estimation is quantified. The medium scenario is selected as the subsequent optimization target. In the missing data sensitivity analysis, the conservative scenario is: missing data = 0, the medium scenario is: missing data = average hourly load, and the high scenario is: missing data = 1.5 times the average hourly load.

5. The method for analyzing the entire chain of watershed water pollution based on emission inventories and water quality monitoring according to claim 1, characterized in that, Step three, specifically the calculation of monitoring load and the establishment of the whole-chain coefficient system, includes the following: (a) Based on hourly monitoring data, the annual pollutant load is calculated using the integral method: ； Where L is the annual pollutant load, kg; C i Q represents the pollutant concentration monitored in the i-th hour, in mg / L; i Let m be the instantaneous flow rate monitored in the i-th hour. 3 / s; i ranges from 1 to n, where n is the number of valid monitoring hours throughout the year, and 3.6 is the unit conversion factor; (b) Establish a three-level coefficient system of "discharge → inflow into the river → monitoring load" and calculate the comprehensive inflow coefficient, river transport coefficient and comprehensive transport coefficient respectively.

6. The method for analyzing the entire chain of watershed water pollution based on emission inventories and water quality monitoring according to claim 5, characterized in that, Comprehensive inflow coefficient = inflow volume / discharge volume; River transport coefficient = monitoring load / inflow volume; Comprehensive transport coefficient = monitoring load / discharge volume.

7. The method for analyzing the entire chain of watershed water pollution based on emission inventories and water quality monitoring according to claim 1, characterized in that, In step four, the specific method for Bayesian global coefficient optimization is as follows: To address the severely underdetermined problem, a Bayesian maximum a posteriori estimation framework is adopted: (a) Problem formalization: ； Among them, M p To supplement the load of the monitoring section, Let p be the emission amount of pollutant from the j-th source. Let p be the initial inflow coefficient of pollutant p from the j-th source. U is the correction factor for the inflow coefficient of pollutant p from the j-th source to be optimized. p Contribute to unknown sources that were not included in the list; j ranges from 1 to m; (b) Differentiated prior distribution: The correction factor uses a truncated normal prior, in the following form: ； Prior mean μ of each source j and variance σ j 2 The settings are differentiated based on data reliability; a and b are the lower and upper limits of the truncated normal distribution, used to constrain the reasonable physical range of the correction factor; Unknown source Using Gamma distribution prior Its shape parameter ν and velocity parameter The design aims to encourage contributions from smaller, unknown sources; (c) Likelihood function: Assuming the observation error follows a normal distribution, the likelihood function is: ； In the formula: For load prediction based on the parameter θ model; σ obs Let M be the standard deviation of the observation error. p 10%; (d) Solution algorithm: The differential evolution global optimization algorithm is used to solve the MAP estimation, and the 95% confidence interval of the parameters is calculated by Laplace approximation; (e) Anomaly source detection: Calculate the correction factor for each source. ,in μ is the optimal estimate of the j-th source correction factor. j and σ j Let z be the mean and standard deviation of its prior distribution; when z j When the threshold is exceeded, the source is marked as an abnormal source, providing a priority ranking for data verification.

8. The method for analyzing the entire chain of watershed water pollution based on emission inventories and water quality monitoring according to claim 1, characterized in that, In step five, the specific method for secondary optimization of functional zones for complex pollution sources is as follows: For urban non-point source pollution sources with complex internal structures, the global correction factor is decomposed to the functional zone level: (a) Model establishment: The formula for calculating urban non-point source pollution inflow into rivers is: ； Among them, W a Cj represents the effective weight of the a-th functional area; p,a denoted as pollutant p in functional zone a, in mg / L; R is the stormwater pipe network correction factor, h is the total number of functional zones, and D is the distance correction factor. The initial values ​​of R and D are determined based on the "Technical Guidelines for Compiling Water Pollutant Source Discharge Inventory" or local experience and are used as parameters to be estimated in the optimization. (b) Parameter optimization: There are a total of 4h+2 parameters and 4 constraint equations. The same Bayesian MAP framework as in step four is used to solve the problem, and the prior distribution center is set as the original coefficient value. (c) Contribution structure redistribution: Analyze the changes in the contribution of each functional zone to the river inflow before and after optimization, and identify the main sources of parameter deviation.

9. The method for analyzing the entire chain of watershed water pollution based on emission inventories and water quality monitoring according to claim 1, characterized in that, Step Six: The specific method for Monte Carlo uncertainty analysis is as follows: (a) Set the emissions of each pollution source to follow a normal distribution and the inflow coefficient to follow a uniform distribution within a multiple range, and conduct N random sampling simulations; (b) Obtain the probability distribution characteristics of the total amount of each pollutant entering the river, and compare it with the monitoring load to assess the degree of closure; (c) Perform elasticity coefficient analysis: Increase the emission amount by a fixed percentage at each source and calculate the proportion of the impact on the total amount entering the river. Identify the most sensitive sources of pollution.

10. The method for analyzing the entire chain of watershed water pollution based on emission inventories and water quality monitoring according to claim 1, characterized in that, The specific methods for establishing the spatial distance-attenuation model and analyzing its effective contribution in step seven are as follows: (a) Spatial decay model: ； Among them, M p To supplement the load of the monitoring section, Let p be the emission amount of pollutant from the j-th source. Let k be the initial inflow coefficient of pollutant p from the j-th source. p The river channel decay rate of pollutant p, in km -1 ;d j γ is the river channel distance from the j-th source to the outflow monitoring section, in km; p is the global correction factor for pollutant p; for area sources, the representative distance is taken by aggregating according to the control unit, and for point sources, the distance to the monitoring section is calculated for each individual source; the value of j is 1~m; (b) Joint optimization: The differential evolution algorithm is used to simultaneously optimize the decay rate k. p and global correction factor γ p This allows the model predictions to approximate the monitoring load after completion; each pollutant is optimized independently. (c) Effective contribution rearrangement: After considering the spatial distance attenuation, the effective contribution rate of each source to the monitoring section is recalculated to reveal the spatial redistribution effect between "emission contribution" and "water quality impact".

Images