Method and system for identifying the contribution of nitrogen and phosphorus sources in a watershed using a data system

Abstract

The present invention provides a method and system for identifying the contribution of nitrogen and phosphorus sources in a watershed using a data system. The method comprises the following steps: Step 1: Collecting water samples during the dry and wet seasons of the target study year for nitrate nitrogen and oxygen isotope determination, and collecting water samples during any season for phosphate oxygen isotope determination, obtaining isotope data at three key time points; Step 2: Estimating the instantaneous source contribution ratio based on the obtained isotope data; Step 3: Converting the estimated instantaneous source contribution ratio into a Dirichlet weak prior distribution; Step 4: Establishing a regression relationship between the source ratio and an annual proxy variable; Step 5: Constructing a four-level Bayesian unified inference model, employing the Markov Chain Monte Carlo method for parameter estimation, and outputting source ratios and loads for multiple years, along with their uncertainties. The present invention utilizes low-cost time-series proxy variables to describe the interannual variability of source ratios, eliminating the long-standing reliance on high-frequency isotope sampling.

Description

Technical Field

[0001] The present invention relates to the technical field of watershed nitrogen and phosphorus source contribution identification, and in particular to a method and system for identifying watershed nitrogen and phosphorus source contributions using a data system. Background Art

[0002] Eutrophication in rivers worldwide is currently increasing to varying degrees. Against this backdrop, reducing the input of exogenous nutrients from land basins into river waters has become a necessary measure for controlling eutrophication. Nitrogen and phosphorus are the two main types of nutrients that contribute to eutrophication. Therefore, controlling the exogenous input of nitrogen and phosphorus nutrients into rivers is a key approach to controlling eutrophication. Currently, concentration indicators for various forms of nitrogen and phosphorus have become a routine component of water environment monitoring in all river basins. However, a single concentration or load indicator is ineffective for precise governance. Identifying the source of excessive concentrations and formulating spatially and temporally differentiated emission reduction targets accordingly are the prerequisites for achieving precise and effective control of exogenous nitrogen and phosphorus nutrient inputs.

[0003] Traditional methods for nitrogen and phosphorus pollution source attribution include emission inventories / empirical coefficients, empirical regression methods (such as SPARROW and WRTDS-Land), process models (such as SWAT, HSPF, and HYPE), and fingerprint tracking / hybrid models. The emission inventory / empirical coefficient method uses statistics on fertilizer, population, and livestock by administrative region, multiplying these by a fixed loss coefficient to directly assign loads to rivers. This method does not require water quality monitoring, but due to the assumption that the coefficients are universal across regions and the method's limitations of ignoring river interception, the spatial and temporal resolution of the pollution source attribution results is coarse, making it unsuitable for dynamic assessment. The empirical regression method uses numerical regression of pollution loads from multiple rivers with statistical values ​​of land use and climate indicators to identify the contributions of individual pollution sources. While this method can omit the investigation of pollution characteristics of polluting end members, it cannot distinguish between fertilizer and soil mineralization, both of which belong to cultivated land. Furthermore, this method is difficult to generalize across regions. The process model method utilizes a hydrological model, dividing the watershed into hydrological response units and simulating the processes of rainfall, runoff, erosion, plant uptake, and denitrification. Although the hydrological-water quality model can realize scenario simulation and measure evaluation, its input parameters are high-dimensional and the data volume is huge, and it has high requirements for monitoring frequency and regional coverage.

[0004] and As naturally occurring stable isotopes in water, different pollutant end-members have their own unique ratio characteristics that serve as their detectable "fingerprints." The isotopic signature of polluted water will appear as a weighted average of the individual pollutant end-members due to mass conservation. Hybrid models such as MixSIAR, simmr, and Dirichlet-LogNormal can utilize statistical methods such as maximum likelihood or Bayesian methods to account for errors and prior characteristics, superimposing different fingerprints into the same hybrid model, improving the discernibility of fingerprint information and enabling quantitative analysis of the contribution ratio of pollution sources. This method has the advantages of accurately distinguishing pollution sources within a single river, being quantitatively interpretable, and being easy to promote across different rivers.

[0005] However, global data on river conduction is sparse 、 The main problems faced by isotope fingerprint / mixture model source apportionment are:

[0006] (1) The monitoring budget is limited and cannot be carried out continuously or frequently 、 Isotope sampling and even water quality monitoring. To monitor the annual load, it is necessary to cover the hydrological fluctuations throughout the year, generally 2 to 4 times a month, or 24 to 48 times a year.

[0007] (2) It is difficult for developing regions to maintain multi-year series, resulting in the inability to obtain information on interannual source structure changes.

[0008] Patent application CN119207601A discloses a method for tracing the source of pollutants in a watershed using multiple technologies and quantifying their load. This method involves pollution monitoring technology and involves sample collection, experimental analysis and data processing, constructing a SWAT model, validating the SWAT model's results, and identifying and quantifying pollution sources using nitrogen and oxygen isotope analysis. However, this patent fails to fully address existing technical issues and does not meet the requirements of the present invention. Summary of the Invention

[0009] In view of the defects in the prior art, the purpose of the present invention is to provide a method and system for identifying the contribution of nitrogen and phosphorus sources in a watershed through a data system.

[0010] The method for identifying the contribution of nitrogen and phosphorus sources in a watershed according to the data system provided by the present invention comprises:

[0011] Step 1: Collect water samples during the dry season of the target study year for nitrate nitrogen and oxygen isotope determination, during the wet season for nitrate nitrogen and oxygen isotope determination, and in any season for phosphate oxygen isotope determination, obtaining isotope data at three key time points in total;

[0012] Step 2: Based on the obtained isotope data, estimate the contribution ratio of the transient source;

[0013] Step 3: Convert the estimated instantaneous source contribution ratio into a Dirichlet weak prior distribution;

[0014] Step 4: Establish a regression relationship between source proportion and annual proxy variables, including cultivated land coverage, nitrogen fertilizer application intensity, night light index, and population density;

[0015] Step 5: Construct a four-level Bayesian unified inference model consisting of an observation layer, a source ratio layer, a regression layer, and a hyperparameter layer, wherein: the observation layer describes the load observation value, the source ratio layer describes the Dirichlet weak prior distribution, the regression layer describes the regression relationship between the source ratio and the annual proxy variable, and the hyperparameter layer defines the prior distribution of the regression layer parameters; the Markov chain Monte Carlo method is used for parameter estimation, and the source ratio and load and their uncertainties for multiple years are output.

[0016] Preferably, the step 1 comprises: the nitrogen isotope determination during the dry season comprises Double isotope determination; the nitrogen isotope determination during the flood season includes Dual isotope determination; the phosphorus isotope determination in any season is Isotope determination;

[0017] The step 2 includes: estimating the instantaneous source contribution ratio by adopting the least square method or the MixSIAR mixed model.

[0018] Preferably, the step 3 includes:

[0019] The estimated instantaneous source contribution ratio Snap is converted to Dirichlet weak prior distribution:

[0020]

[0021] in, is the precision parameter; is a weak prior distribution function;

[0022] according to Construct Dirichlet distribution. The specific construction steps are as follows:

[0023] Calculate the scaled source ratio :

[0024]

[0025] in, K is the number of pollution sources, is through the precision parameter The adjusted proportion of each source, ;

[0026] all After scaling, take the absolute value to get the weight of each source, and then build the model through Dirichlet distribution:

[0027]

[0028] The final Dirichlet prior distribution is used for subsequent Bayesian inference. In the Bayesian framework, this distribution will be combined with subsequent data through Bayes' theorem to derive the posterior distribution of the source proportion.

[0029] Preferably, step 4 includes:

[0030] Establish a regression relationship between source proportion and proxy variable:

[0031]

[0032] in, is the logistic transformation of the source ratio, which represents the logarithmic ratio of the source ratio; is the intercept term of the regression model, which represents the estimated value of the dependent variable when the independent variable is zero; is the parameter to be estimated; is the set of independent variables; is the error term.

[0033] Preferably, the step 5 includes:

[0034] Describe the relationship between observed data and model parameters:

[0035]

[0036] in, is the water quality observation value, is the design matrix, is the parameter to be estimated, is the variance of the error; is the normal distribution function;

[0037] Use a Dirichlet distribution for the source ratio:

[0038]

[0039] Select initial values ​​for the various model parameters. The source scale parameter is initially set to a uniform distribution:

[0040]

[0041] In each iteration, the parameters are updated by the Metropolis-Hastings algorithm or Gibbs sampling;

[0042] After the run, extract the posterior distribution of each parameter and use the posterior mean To estimate the relative contribution of the sources:

[0043]

[0044] in, is the actual water quality observation dataset, is the expected value of the posterior distribution.

[0045] The system for identifying the contribution of nitrogen and phosphorus sources in a watershed according to the data system provided by the present invention includes:

[0046] Module M1: Collect water samples for nitrate nitrogen and oxygen isotope determination during the dry season of the target study year, for nitrate nitrogen and oxygen isotope determination during the wet season, and for phosphate oxygen isotope determination in any season, obtaining isotope data at three key time points.

[0047] Module M2: Estimate the contribution ratio of transient sources based on the obtained isotope data;

[0048] Module M3: Convert the estimated instantaneous source contribution ratio into a Dirichlet weak prior distribution;

[0049] Module M4: Establishing a regression relationship between source ratio and annual proxy variables, including cultivated land coverage, nitrogen fertilizer application intensity, night light index, and population density;

[0050] Module M5: Construct a four-level Bayesian unified inference model including observation layer, source ratio layer, regression layer and hyperparameter layer, wherein: the observation layer describes the load observation value, the source ratio layer describes the Dirichlet weak prior distribution, the regression layer describes the regression relationship between the source ratio and the annual proxy variable, and the hyperparameter layer defines the prior distribution of the regression layer parameters; the Markov chain Monte Carlo method is used for parameter estimation, and the source ratio and load of multiple years and their uncertainties are output.

[0051] Preferably, the module M1 includes: the nitrogen isotope determination during the dry season includes Double isotope determination; the nitrogen isotope determination during the flood season includes Dual isotope determination; the phosphorus isotope determination in any season is Isotope determination;

[0052] The module M2 includes: estimating the instantaneous source contribution ratio by adopting the least square method or the MixSIAR mixed model.

[0053] Preferably, the module M3 includes:

[0054] The estimated instantaneous source contribution ratio Snap is converted to Dirichlet weak prior distribution:

[0055]

[0056] in, is the precision parameter; is a weak prior distribution function;

[0057] according to Construct Dirichlet distribution. The specific construction steps are as follows:

[0058] Calculate the scaled source ratio :

[0059]

[0060] in, K is the number of pollution sources, is through the precision parameter The adjusted proportion of each source, ;

[0061] all After scaling, take the absolute value to get the weight of each source, and then build the model through Dirichlet distribution:

[0062]

[0063] The final Dirichlet prior distribution is used for subsequent Bayesian inference. In the Bayesian framework, this distribution will be combined with subsequent data through Bayes' theorem to derive the posterior distribution of the source proportion.

[0064] Preferably, the module M4 includes:

[0065] Establish a regression relationship between source proportion and proxy variable:

[0066]

[0067] in, is the logistic transformation of the source ratio, which represents the logarithmic ratio of the source ratio; is the intercept term of the regression model, which represents the estimated value of the dependent variable when the independent variable is zero; is the parameter to be estimated; is the set of independent variables; is the error term.

[0068] Preferably, the module M5 includes:

[0069] Describe the relationship between observed data and model parameters:

[0070]

[0071] in, is the water quality observation value, is the design matrix, is the parameter to be estimated, is the variance of the error; is the normal distribution function;

[0072] Use a Dirichlet distribution for the source ratio:

[0073]

[0074] Select initial values ​​for the various model parameters. The source scale parameter is initially set to a uniform distribution:

[0075]

[0076] In each iteration, the parameters are updated by the Metropolis-Hastings algorithm or Gibbs sampling;

[0077] After the run, the posterior distribution of each parameter is extracted and the posterior mean is used to estimate the relative contribution of the source:

[0078]

[0079] in, is the actual water quality observation dataset, is the expected value of the posterior distribution.

[0080] Compared with the prior art, the present invention has the following beneficial effects:

[0081] (1) Reduce monitoring costs and improve source attribution accuracy. By constructing a Dirichlet weak prior distribution using isotope snapshots at key time points (nitrogen isotopes during the dry / wet seasons + phosphorus and oxygen isotopes during any season), the reasonable range of pollution source contribution ratios can be limited (physical anchoring). Low-cost time series proxy variables (cultivated land coverage, nitrogen fertilizer intensity, nighttime lights, and population density) can be used to describe the interannual variation in source proportions, thereby breaking away from the long-term reliance on high-frequency isotope sampling and improving source attribution accuracy.

[0082] (2) Completely propagate uncertainty and avoid overfitting. Construct a four-layer Bayesian model (observation layer, source ratio layer, regression layer, and hyperparameter layer), and use the MCMC algorithm to uniformly infer multi-year source ratios and loads to fully propagate uncertainty; introduce the Dirichlet weak prior (accuracy parameter ), balances prior constraints and data adaptability, avoids overfitting, and improves the generalization ability of the model.

[0083] (3) Support long-term trend analysis and easy to promote and apply. Integrate global public time series data (ESA-CCI cultivated land, VIIRS night lights, WorldPop population, etc.), replace traditional high-cost continuous monitoring, and support long-term trend analysis. BRIEF DESCRIPTION OF THE DRAWINGS

[0084] Other features, objects and advantages of the present invention will become more apparent upon reading the detailed description of non-limiting embodiments with reference to the following drawings:

[0085] Figure 1 Flowchart of the method for identifying the contribution of nitrogen and phosphorus sources in a watershed for the data system of the present invention. DETAILED DESCRIPTION

[0086] The present invention will be described in detail below with reference to specific embodiments. The following examples will help those skilled in the art to further understand the present invention, but are not intended to limit the present invention in any form. It should be noted that, for those skilled in the art, several changes and improvements can be made without departing from the scope of the present invention. These all fall within the scope of protection of the present invention.

[0087] Example 1

[0088] like Figure 1 The present invention provides a method for identifying the contribution of nitrogen and phosphorus sources in a watershed using a data system, comprising:

[0089] 1. Snapshot Sampling Strategy

[0090] Implement a streamlined sampling plan during the target study year, including:

[0091] (1) Dry season Dual isotope determination;

[0092] (2) Flood season Dual isotope determination;

[0093] (3) Any season Isotope determination;

[0094] Based on the isotope data at the three key time points mentioned above, the contribution ratio of the transient source was estimated using the least squares method or the MixSIAR mixing model. snap.

[0095] 1. The specific process of isotope determination is as follows:

[0096] (1) Sample collection: Based on the above design, select appropriate time periods (dry and wet seasons) and collect water samples from representative locations. After collecting the samples, store them in appropriate containers and analyze them quickly to minimize sample changes.

[0097] (2) Sample pretreatment:

[0098] Nitrogen extraction: Ion exchange method was used to extract nitrogen from water samples. , remove impurities and concentrate the sample.

[0099] Phosphorus extraction: similar operation is used to extract from water samples .

[0100] (3) Isotope analysis: gas chromatography-mass spectrometry (GC-IRMS) Determination and This method can simultaneously analyze isotope ratios.

[0101] Using isotope mass spectrometry Isotope analysis was performed to determine the value.

[0102] (4) Data processing: Record the isotope ratio of each sample to ensure the accuracy and traceability of the data.

[0103] (5) Result recording: Organize and archive the measured isotope ratios to prepare for subsequent analysis.

[0104] 2. Expression and estimation process of least squares method or MixSIAR mixture model;

[0105] The expression of the least squares method is:

[0106]

[0107] in: n is the total number of samples; It is i The observed value of a sample; is the independent variable value corresponding to each sample; are the parameters to be estimated.

[0108] The source ratio estimation expression of the MixSIAR mixing model is:

[0109]

[0110] in: k is the number of pollution sources; It is j Contribution ratio of individual sources; It isi In the sample j Isotope values ​​of individual sources; is the error term.

[0111] Estimation process:

[0112] 1. Data preparation: Collect isotope data and related indicators corresponding to water samples and establish a data set;

[0113] 2. Model construction: Build a parameter model based on the MixSIAR hybrid model and set the objective function;

[0114] 3. Run the model: solve the optimal source ratio through numerical optimization method i ;

[0115] 4. Result verification: Test the model’s fitting effect to ensure that it can accurately reflect the source proportion.

[0116] 2. Construction of weak prior distribution

[0117] Convert the source proportion estimate obtained by snapshot sampling into a Dirichlet weak prior distribution:

[0118]

[0119] The accuracy parameter Setting it to 0.03–0.05 corresponds to a statistical information size of approximately 2–3 valid observations. This parameter setting ensures that the prior information provides reasonable constraints rather than excessive restrictions, allowing the posterior distribution to dynamically adjust based on the observed data.

[0120] Overview of steps:

[0121] 1. Calculate source ratio

[0122] The isotope data of water samples were analyzed using the least squares method or MixSIAR mixing model to obtain the estimated proportion of each pollution source. snap .

[0123] Specifically, snap By analyzing the sampling data, we can determine the percentage of each pollution source that contributes to the nitrogen, phosphorus and other components in the water sample. For example, if there are three main pollution sources A, B, and C, we may get , indicating the contribution ratio of sources A, B and C respectively.

[0124] 2. Set the precision parameters ( )

[0125] Choose the appropriate precision parameters , this parameter is usually set between 0.03 and 0.05.

[0126] This parameter reflects the degree of trust in prior information. The smaller the value, the less constraint on the prior and the greater flexibility the model allows; while the larger the value, the stronger the dependence on existing information.

[0127] 3. Constructing Dirichlet Distribution

[0128] According to the adjusted ratio To construct the Dirichlet distribution. The Dirichlet distribution is a multidimensional probability distribution that is suitable for representing the probability distribution of a random vector consisting of multiple positive numbers whose sum is 1.

[0129] The specific construction steps are as follows:

[0130] Calculate the scaled source ratio:

[0131]

[0132] in K is the number of pollution sources, is through the precision parameter The adjusted proportion of each source, i =1,2,…, K .

[0133] Make sure all After scaling correctly, take the absolute value (because snap The components of must be greater than 0), and the weight of each source is obtained.

[0134] The model can then be constructed using the Dirichlet distribution:

[0135] Dirichlet , ,…, )

[0136] The resulting Dirichlet prior distribution is used in subsequent Bayesian inference. In the Bayesian framework, this distribution is combined with subsequent data using Bayes' theorem to derive the posterior distribution of source proportions, enabling more accurate pollution source analysis.

[0137] 3. Modeling of spatiotemporal variability

[0138] Taking into account the interannual changes in the pollution source structure, a regression relationship between the source ratio and the proxy variable is established:

[0139]

[0140] in, It is the logistic transformation of the source proportion, which represents the logarithmic ratio of the source proportion and is often used for modeling in regression analysis; is the intercept term of the regression model, which represents the estimated value of the dependent variable when the independent variable is zero; It is the coefficient in the regression model, which represents the expected change in the dependent variable when the independent variable changes by one unit; is a set of independent variables, including indicators such as cultivated land coverage, nitrogen fertilizer application intensity, night light index and population density; is the error term, which represents the unexplained part of the model or random fluctuations.

[0141] in include:

[0142] (1) Cultivated land coverage (based on remote sensing data);

[0143] (2) Nitrogen fertilizer application intensity (NAPPN dataset);

[0144] (3) Night light index (VIIRS data);

[0145] (4) Population density (WorldPop dataset);

[0146] These global public datasets provide continuous time series from 1992 to the present, supporting long-term trend analysis and forecasting.

[0147] 4. Hierarchical Bayesian Inference Framework

[0148] Construct a four-level Bayesian unified inference model:

[0149] Level 1 - Observation layer: load observations Follows a lognormal distribution;

[0150] Level 2 - Source Scale Layer: ;

[0151] Level 3 - Regression Layer: ;

[0152] Level 4 - Hyperparameters layer: ;

[0153] The Markov Chain Monte Carlo (MCMC) method is used for parameter estimation, and the source ratio for the period of 5 to 10 years is also output. and load , and fully propagate uncertainty.

[0154] MCMC is a powerful statistical tool, particularly suitable for estimating parameters of complex models in Bayesian analysis. The specific steps are as follows:

[0155] 1. Model Construction

[0156] Building a complete Bayesian model: Based on the research question, construct a Bayesian model that includes an observation layer, a parameter layer, and a prior distribution layer. Observation layer: describes the relationship between the observed data and the model parameters. For example, suppose we have water quality observations , we can express it as follows:

[0157]

[0158] in, is the design matrix, is the parameter to be estimated, is the variance of the error.

[0159] Set prior distribution: Introduce prior knowledge into the model. Common prior distributions include normal distribution, Gamma distribution, Dirichlet distribution, etc. For example, Dirichlet distribution can be used for source ratio:

[0160]

[0161] 2. Chain initialization

[0162] Select initial values: Choose reasonable initial values ​​for each model parameter. When choosing initial values, ensure that they are reasonable to avoid the influence of the initial state of the chain on the estimation results. For example, the source scale parameter can be initially set to a uniform distribution:

[0163]

[0164] Start multiple chains: It is recommended to start at least three independent chains to improve the robustness of the results. Each chain starts with a different initial value, which helps avoid local optimal solutions.

[0165] 3. Iterative Updates

[0166] Parameter update: In each iteration, the parameters are updated by the Metropolis-Hastings algorithm or Gibbs sampling.

[0167] Metropolis-Hastings algorithm: randomly generates new parameter values ​​and calculates their posterior probabilities. The acceptance rate determines whether to accept the new value:

[0168] (1) Generate a new parameter value , typically generated randomly within a small neighborhood of the current location.

[0169] (2) Calculate the current parameter value and new parameter values The posterior probability of :

[0170]

[0171] (3) Decide whether to accept based on the acceptance rate α:

[0172]

[0173] Gibbs sampling: Sample each parameter in turn. Sampling is performed under the condition that other parameters are fixed. For example, if there are two parameters in the model and , then alternate sampling:

[0174] (1) From the conditional distribution Sampling new .

[0175] (2) Then from Sampling new .

[0176] 4. Convergence check

[0177] Observe chain trajectories: During the operation of multiple chains, periodically check the chain trajectory graph. Merge the samples of all chains. If they fluctuate on the same plane and cover similar areas, it means that the chains are gradually converging.

[0178] Calculate R-hat: Using the Gelman-Rubin diagnostic, an R-hat value less than 1.1 generally indicates good convergence of the model. The calculation method is as follows:

[0179] (1) Calculate the sample mean and variance of each chain.

[0180] (2) Combine all chain results to get the overall mean and variance.

[0181] (3) Calculate the degree of dispersion and compare it with the internal differentiation of each chain to obtain R-hat.

[0182] Valid sample number: Check the valid sample number of each parameter, which is usually required to be greater than 1000 to ensure the stability and accuracy of the estimation.

[0183] 5. Output results

[0184] Generate posterior distributions: After the run, extract the posterior distributions for each parameter. These posterior distributions represent the model's latest understanding of source proportions and load estimates. For example, the posterior means can be used to estimate the relative contributions of the sources:

[0185] .

[0186] Example 2

[0187] Example 2 is a preferred example of Example 1.

[0188] The present invention provides a method for identifying the contribution of nitrogen and phosphorus sources in a watershed using a data system, comprising:

[0189] 1. Sampling Design

[0190] a. Take one sample during the dry season: ;

[0191] b. Take one sample during the flood season: ;

[0192] c. Take one sample from any season: ;

[0193] A total of 3 samples completed the snapshot for the year.

[0194] 2. Snapshot Endmember Solution

[0195] (1) Collect or measure end-member isotopes μ±SD (fertilizer, sewage, soil, etc.).

[0196] (2) Use MixSIAR / least squares to find snap , and record the covariance Σ.

[0197] 3. Constructing Dirichlet Weak Prior

[0198] (1) Assume κ = 0.04 (which can be automatically adjusted according to the sample size).

[0199] (2) Calculation , as the Level-2 prior parameter in the hierarchical model.

[0200] 4. Organizing annual driving variables

[0201] (1) During the study period, the following data were obtained each year: ESA-CCI cropland%, NAPPN remote sensing nitrogen fertilizer, VIIRS night light, WorldPop population, etc.

[0202] (2) Import form And standardized to mean 0 and variance 1.

[0203] 5. Load Estimation

[0204] Use monthly water quality and flow rates to run WRTDS or loadReg to get an annual load estimate and error .

[0205] The main process includes:

[0206] 1. WRTDS model construction and load estimation

[0207] (1) Model construction:

[0208] The WRTDS model is a nonparametric regression method used to establish the relationship between concentration (C) and flow (Q), time (t) and season (season). The basic form of the WRTDS model is:

[0209]

[0210] in: is the concentration; For traffic; is the time variable; Used to express seasonal effects in models; Denotes the error term. Determines the model window width (such as time, flow, and seasonal windows). These parameters will affect the smoothness and flexibility of the model.

[0211] (2) Model fitting:

[0212] The WRTDS model was fitted using the collected monthly water quality and flow data to obtain predicted concentrations.

[0213] 2. Load calculation

[0214] (1) Monthly load estimation: Based on the concentration predicted by the WRTDS model and the actual flow data, the monthly load is calculated using the following formula:

[0215]

[0216] in: is the load in month t (e.g. kg / month); is the predicted concentration for that month (in mg / L); is the flow rate of the month (in m³ / s); k is the unit conversion factor used to ensure consistent load units.

[0217] (2) Annual load estimation: Sum the loads of all months to obtain the annual load estimate:

[0218]

[0219] 3. Uncertainty estimation ( )

[0220] (1) Uncertainty estimation methods:

[0221] Bootstrap resampling:

[0222] Randomly sample the original data and repeat this process to generate multiple resampled datasets of equal size. The model is fitted to each resampled dataset, and the annual load is calculated. After multiple resamplings, the standard deviation of the annual loads obtained from all resamples is calculated as an indicator of uncertainty.

[0223] (2) Error formula:

[0224] The error in load estimation can be expressed as:

[0225]

[0226] in: It is i The actual load of the observation; is the load predicted by the WRTDS model.

[0227] 6. Assembling the Stan / JAGS model

[0228] 1. Model Construction

[0229] Use Bayesian methods such as Stan or JAGS to construct a Dirichlet-LogNormal model. This model typically consists of four layers: the observation layer, the source ratio layer, the regression layer, and the hyperparameter layer. This hierarchical structure allows the model to effectively incorporate observational data, source ratio information, and individual effects.

[0230] (1) Observation layer: describes the observation load ( ) and the model predicted value. Normal distribution is often used to represent the observed load. The basic form of this layer is:

[0231]

[0232] Where μ is the predicted value of the load and σ is the standard deviation of the error.

[0233] (2) Source ratio layer: Dirichlet distribution is used to model the proportional relationship between different pollution sources and the total load. Its expression is:

[0234]

[0235] Where θ is the source ratio vector, K is the number of pollution sources, is a hyperparameter associated with the jth pollution source.

[0236] (3) Regression layer: The regression layer expresses the impact of independent variables (such as cultivated land coverage, fertilization intensity, etc.) on load. The linear regression model can be used to express the regression relationship:

[0237]

[0238] Where X is the independent variable, β is the regression coefficient, and α is the intercept.

[0239] (4) Hyperparameter layer:

[0240] Hyperparameters are used to define prior information and allow updating parameters in the model. This layer is an important part of ensuring the reasonableness of the source ratio distribution and error distribution.

[0241] 2. Use Stan model construction

[0242] The modeling process in Stan is as follows. The modeling process in JAGS is similar to that in Stan:

[0243] (1) Data block: defines the data required for the model, including observation load, independent variables and sample size. For example: N represents the number of samples; represents the observed load; X Represents the independent variable.

[0244] (2) Parameter block: defines the parameters to be estimated, such as source ratio i and standard deviation of error s .

[0245] (3) Model block: Define the relationship in the model block. First, describe the load observation model, and then set the Dirichlet prior of the source ratio, which is usually assumed to be:

[0246]

[0247] in, kappa is a hyperparameter that controls the strength of the prior.

[0248] 3. Model calculation process

[0249] The model operation process mainly includes the following steps:

[0250] (1) Model implementation: The observation load in the data set is converted into , independent variables X Verify that the data is in the correct format and size to meet the model requirements.

[0251] (2) Model fitting: Run the MCMC process to update the model parameters. MCMC methods can effectively sample in a higher-dimensional parameter space, thereby obtaining parameter estimates from the posterior distribution.

[0252] (3) Result analysis: Obtain the posterior distribution of the model parameters. Parameters can be analyzed using statistical measures (such as mean, variance, and confidence interval). For example, a 95% confidence interval can be used to assess the credibility of the parameters.

[0253] 7. Result Extraction and Diagnosis

[0254] (1) Output posterior mean and 95% CI: .

[0255] (2) Check that R-hat < 1.1 and the number of valid samples > 1000, and draw a trace-plot.

[0256] 8. Error Assessment (Optional Cross-Validation)

[0257] Error assessment is an important part of model training and verification. Its purpose is to evaluate the accuracy and stability of the model through a series of quantitative indicators.

[0258] 1. Extrapolation Verification Method

[0259] Retain one or rolling year, extrapolation validation is to test the model's predictive ability for unknown data by retaining part of the data set. The specific method is as follows:

[0260] (1) Data division:

[0261] Divide the entire dataset into two parts: a training set and a test set. You can choose to divide the data by year, retaining one year's data for the test set and the remaining years' data for the training set. Another approach is to use a rolling window strategy, retaining a portion of the data for validation over multiple years.

[0262] (2) Model training: Use the training set data to fit the model and obtain the model parameters. The data in the training set should generally contain enough time points and samples to ensure that the model can learn the underlying trends and relationships in the data.

[0263] (3) Use the fitted model to predict the test set data and calculate the model's load prediction value. Record the actual observed load value of each retained sample and the model prediction value for subsequent calculation of nRMSE and coverage.

[0264]

[0265] Where X is the independent variable, β is the regression coefficient, and α is the intercept.

[0266] (4) Performance index calculation:

[0267] Normalized Root Mean Square Error (nRMSE): Calculating nRMSE is a standard way to measure the accuracy of a model’s predictions and is defined as:

[0268]

[0269] in, represents the predicted value of the i-th sample, represents the corresponding observation value. The closer the value is to 0, the better the predictive ability of the model.

[0270] Coverage: Calculate the confidence intervals for the load forecast (usually 95%) and check whether these confidence intervals cover the actual observed values. The specific method is to record the proportion of the actual observed load in the test set that falls within the confidence interval, which is defined as coverage. The calculation formula is:

[0271]

[0272] Where N is the number of test samples, Is an indicator function that indicates whether the observed value falls within the predicted confidence interval .

[0273] 2. Error analysis and adjustment

[0274] If the error is greater than the preset threshold, after extrapolation verification, if the value of nRMSE or coverage exceeds the reasonable threshold set by the research institute, the following steps need to be taken to adjust the model:

[0275] (1) Adjustment value: The value is the prior strength parameter of the source ratio, which controls the shape of the Dirichlet distribution. If the nRMSE value is too large, it may mean that the model is not accurate enough in estimating the source ratio. In this case, you can reduce the Adjusting the value of α increases the model's sensitivity to the existing data, causing the model to "follow" the training data more closely. This typically results in greater variance and less bias, leading to improved predictive power.

[0276] (2) Introducing high-weight proxy variables: If the adjustment If the model's predictive performance still fails to improve after adding the variables, consider introducing high-weighted proxy variables. For example, soil moisture, precipitation, and fertilizer application rates can affect nitrogen and phosphorus losses and concentrations. The impact of high-weighted proxy variables can be confirmed through preliminary analysis and domain knowledge to ensure their importance in the analysis.

[0277] (3) Retraining the model: After adjusting the parameters, return to the model training stage and train the model again using the new parameters and possible new data.

[0278] (4) Continuous monitoring and feedback: After the model is adjusted, continuous monitoring and verification should be carried out to ensure that the effect of each model improvement can be verified by new observation data. Cross-year or cross-regional verification methods can be used to ensure the robustness of the model.

[0279] Those skilled in the art will appreciate that, in addition to implementing the system, device, and various modules provided by the present invention in purely computer-readable program code, it is entirely possible to implement the same program in the form of logic gates, switches, application-specific integrated circuits, programmable logic controllers, embedded microcontrollers, and the like by logically programming the method steps. Therefore, the system, device, and various modules provided by the present invention can be considered a hardware component, and the modules included therein for implementing various programs can also be considered structures within the hardware component; the modules for implementing various functions can also be considered both software programs for implementing the method and structures within the hardware component.

[0280] The above describes specific embodiments of the present invention. It should be understood that the present invention is not limited to the specific embodiments described above, and those skilled in the art may make various changes or modifications within the scope of the claims, which do not affect the essence of the present invention. The embodiments of this application and the features in the embodiments may be combined with each other in any manner unless there is a conflict.

Claims

1. A method for identifying the contribution of nitrogen and phosphorus sources in a watershed using a data system, characterized in that: include: Step 1: Collect water samples during the dry season of the target study year for nitrate nitrogen and oxygen isotope determination, during the wet season for nitrate nitrogen and oxygen isotope determination, and in any season for phosphate oxygen isotope determination, obtaining isotope data at three key time points in total; Step 2: Based on the obtained isotope data, estimate the contribution ratio of the transient source; Step 3: Convert the estimated instantaneous source contribution ratio into a Dirichlet weak prior distribution; Step 4: Establish a regression relationship between source proportion and annual proxy variables, including cultivated land coverage, nitrogen fertilizer application intensity, night light index, and population density; Step 5: Construct a four-level Bayesian unified inference model consisting of an observation layer, a source ratio layer, a regression layer, and a hyperparameter layer, wherein: the observation layer describes the load observation value, the source ratio layer describes the Dirichlet weak prior distribution, the regression layer describes the regression relationship between the source ratio and the annual proxy variable, and the hyperparameter layer defines the prior distribution of the regression layer parameters; the Markov chain Monte Carlo method is used for parameter estimation, and the source ratio and load and their uncertainties for multiple years are output.

2. The method for identifying the contribution of nitrogen and phosphorus sources in a watershed using a data system according to claim 1, characterized in that: The step 1 includes: the nitrogen isotope determination during the dry season includes Double isotope determination; the nitrogen isotope determination during the flood season includes Dual isotope determination; the phosphorus isotope determination in any season is Isotope determination; The step 2 includes: estimating the instantaneous source contribution ratio by adopting the least square method or the MixSIAR mixed model.

3. The method for identifying the contribution of nitrogen and phosphorus sources in a watershed using a data system according to claim 1, characterized in that: The step 3 includes: The estimated instantaneous source contribution ratio Snap is converted to Dirichlet weak prior distribution: in, is the precision parameter; is a weak prior distribution function; according to Construct Dirichlet distribution. The specific construction steps are as follows: Calculate the scaled source ratio : in, K is the number of pollution sources, is through the precision parameter The adjusted proportion of each source, ; all After scaling, take the absolute value to get the weight of each source, and then build the model through Dirichlet distribution: The final Dirichlet prior distribution is used for subsequent Bayesian inference. In the Bayesian framework, this distribution will be combined with subsequent data through Bayes' theorem to derive the posterior distribution of the source proportion.

4. The method for identifying the contribution of nitrogen and phosphorus sources in a watershed using a data system according to claim 1, wherein: The step 4 comprises: Establish a regression relationship between source proportion and proxy variable: in, is the logistic transformation of the source ratio, which represents the logarithmic ratio of the source ratio; is the intercept term of the regression model, which represents the estimated value of the dependent variable when the independent variable is zero; is the parameter to be estimated; is the set of independent variables; is the error term.

5. The method for identifying the contribution of nitrogen and phosphorus sources in a watershed using a data system according to claim 1, characterized in that: The step 5 comprises: Describe the relationship between observed data and model parameters: in, is the water quality observation value, is the design matrix, is the parameter to be estimated, is the variance of the error; is the normal distribution function; Use a Dirichlet distribution for the source ratio: Select initial values ​​for the various model parameters. The source scale parameter is initially set to a uniform distribution: In each iteration, the parameters are updated by the Metropolis-Hastings algorithm or Gibbs sampling; After the run, extract the posterior distribution of each parameter and use the posterior mean To estimate the relative contribution of the sources: in, is the actual water quality observation dataset, is the expected value of the posterior distribution.

6. A system for identifying the contribution of nitrogen and phosphorus sources in a watershed using a data system, characterized in that: include: Module M1: Collect water samples for nitrate nitrogen and oxygen isotope determination during the dry season of the target study year, for nitrate nitrogen and oxygen isotope determination during the wet season, and for phosphate oxygen isotope determination in any season, obtaining isotope data at three key time points. Module M2: Estimate the contribution ratio of transient sources based on the obtained isotope data; Module M3: Convert the estimated instantaneous source contribution ratio into a Dirichlet weak prior distribution; Module M4: Establishing a regression relationship between source ratio and annual proxy variables, including cultivated land coverage, nitrogen fertilizer application intensity, night light index, and population density; Module M5: Construct a four-level Bayesian unified inference model including observation layer, source ratio layer, regression layer and hyperparameter layer, wherein: the observation layer describes the load observation value, the source ratio layer describes the Dirichlet weak prior distribution, the regression layer describes the regression relationship between the source ratio and the annual proxy variable, and the hyperparameter layer defines the prior distribution of the regression layer parameters; the Markov chain Monte Carlo method is used for parameter estimation, and the source ratio and load of multiple years and their uncertainties are output.

7. The system for identifying the contribution of nitrogen and phosphorus sources in a watershed according to claim 6, characterized in that: The module M1 includes: the nitrogen isotope determination during the dry season includes Double isotope determination; the nitrogen isotope determination during the flood season includes Dual isotope determination; the phosphorus isotope determination in any season is Isotope determination; The module M2 includes: estimating the instantaneous source contribution ratio by adopting the least square method or the MixSIAR mixed model.

8. The system for identifying the contribution of nitrogen and phosphorus sources in a watershed according to claim 6, characterized in that: The module M3 includes: The estimated instantaneous source contribution ratio Snap is converted to Dirichlet weak prior distribution: in, is the precision parameter; is a weak prior distribution function; according to Construct Dirichlet distribution. The specific construction steps are as follows: Calculate the scaled source ratio : in, K is the number of pollution sources, is through the precision parameter The adjusted proportion of each source, ; all After scaling, take the absolute value to get the weight of each source, and then build the model through Dirichlet distribution: The final Dirichlet prior distribution is used for subsequent Bayesian inference. In the Bayesian framework, this distribution will be combined with subsequent data through Bayes' theorem to derive the posterior distribution of the source proportion.

9. The system for identifying the contribution of nitrogen and phosphorus sources in a watershed according to claim 6, characterized in that: The module M4 includes: Establish a regression relationship between source proportion and proxy variable: in, is the logistic transformation of the source ratio, which represents the logarithmic ratio of the source ratio; is the intercept term of the regression model, which represents the estimated value of the dependent variable when the independent variable is zero; is the parameter to be estimated; is the set of independent variables; is the error term.

10. The system for identifying the contribution of nitrogen and phosphorus sources in a watershed according to claim 6, characterized in that: The module M5 includes: Describe the relationship between observed data and model parameters: in, is the water quality observation value, is the design matrix, is the parameter to be estimated, is the variance of the error; is the normal distribution function; Use a Dirichlet distribution for the source ratio: Select initial values ​​for the various model parameters. The source scale parameter is initially set to a uniform distribution: In each iteration, the parameters are updated by the Metropolis-Hastings algorithm or Gibbs sampling; After the run, the posterior distribution of each parameter is extracted and the posterior mean is used to estimate the relative contribution of the source: in, is the actual water quality observation dataset, is the expected value of the posterior distribution.

Images