Phosphorus pollution fine source analysis method and system, and storage medium

Abstract

The application provides a phosphorus pollution fine source analysis method and system and a storage medium to solve the problem of inaccurate source analysis caused by the over-reliance of the existing receptor model on the tracer conservation hypothesis and the lack of quantitative characterization means for phosphorus form transformation mechanism, and belongs to the technical field of environmental engineering and water environment monitoring. The method comprises the following steps: collecting phosphorus form abundance data, driving factor data and microbial characteristic data in samples of each receptor sampling point of the river and potential pollution sources; calculating the macro source-sink effect coefficient of each phosphorus form; calculating the biochemical effect coefficient of each phosphorus form; selecting an effective tracer and extracting the transformed abundance parameter of the effective tracer; reconstructing a dynamic correction source spectrum matrix; outputting the final contribution rate distribution of each pollution source to a specific phosphorus form in the receptor; and performing standardization processing and significance test on the basic residual error to identify significant abnormal points.

Description

Technical Field

[0001] This invention belongs to the field of environmental engineering and water environment monitoring technology, specifically a method, system and storage medium for refined source apportionment of phosphorus pollution. Background Technology

[0002] Phosphorus is a key controlling factor in surface water eutrophication, and accurate source apportionment of phosphorus pollution in watersheds is a scientific prerequisite for formulating targeted emission reduction strategies and precise pollution control. Currently, phosphorus pollution source apportionment mainly relies on a combination of stable isotope techniques and receptor models. However, since phosphorus exists in nature with only one available stable isotope, phosphate oxygen, its tracer information capacity is limited. In watershed environments with highly heterogeneous land use and complex pollution source compositions, single isotope fingerprints often face the problem of overlapping characteristic signals, making it difficult to distinguish different pollution sources with similar contribution rates. This leads to significant technical bottlenecks in achieving high spatial resolution and high-precision apportionment using existing methods.

[0003] Traditional phosphorus source apportionment models, such as the Positive Definite Matrix Factorization (PMF) model and the Bayesian Mixture Analysis (MixSIAR) model, largely rely on the conservative assumption of tracers or empirically compensate for deviations in the migration process by introducing static fractionation coefficients. For example, Chinese patent publications CN121141795A and CN121393591A disclose source apportionment methods for river phosphorus pollution based on receptor models, but both assume that phosphorus isotopic fractionation and speciation follow a constant ratio or remain relatively physicochemically stable during its emission from the pollution source to the receptor. However, in real-world complex aquatic environments, phosphorus speciation exhibits highly dynamic and non-conservative behavior, and its migration and transformation are easily affected by the superposition of multiple factors, including complexation of dissolved organic matter, microbial assimilation and decomposition, and fluctuations in local environmental gradients. Although some receptor model frameworks allow for the input of fractionation parameters for correction, existing technologies lack quantitative characterization methods for these complex transformation mechanisms. Because the model relies solely on empirical constants or static prior parameters, it cannot accurately capture the dynamic transformation effects that change with spatial transport and biological disturbance. This makes it easy for the measured signal at the receptor end to deviate from the statically corrected prior source spectrum range when dealing with highly reactive phosphorus pollution scenarios. Consequently, it causes systematic bias in the source apportionment results, resulting in the model failing to converge or outputting posterior probabilities that have no physical meaning.

[0004] Furthermore, most existing analytical methods focus on total phosphorus or single soluble orthophosphate components, failing to effectively extract and utilize multidimensional characterization information of phosphorus speciation. Phosphorus speciation not only reflects the material composition of the initial emission source but also records the physicochemical trajectory of pollutants during migration and transformation. Notably, the non-conservative transformation of phosphorus speciation in complex aquatic environments does not occur in isolation but is highly correlated with the biochemical cycle of dissolved organic matter (DOM). DOM can directly alter the physicochemical state of phosphorus through complexation and competitive adsorption; its molecular components with different degradation characteristics can significantly drive the metabolic activity of the underlying microbial community, thereby dominating the biogeochemical reconstruction process of phosphorus speciation. Given that DOM substantially controls phosphorus speciation transformation at the microscopic level, incorporating it as the core of quantifying non-conservative behavior into analytical models is crucial. Therefore, constructing a refined model with dynamically corrected source spectra by coupling the multidimensional distribution characteristics of phosphorus speciation with the DOM-mediated transformation mechanism is a necessary approach to address the uncertainty of phosphorus source attribution and improve analytical accuracy in current complex watershed environments. Therefore, in response to current technological bottlenecks, developing a source apportionment method that can quantify non-conservative processes and deeply integrate multidimensional characteristics of phosphorus speciation is of great practical value and technical guiding significance for achieving refined attribution of pollution sources in the phosphorus speciation dimension and establishing a precise prevention and control system for phosphorus pollution at the watershed scale. Summary of the Invention

[0005] The technical problem to be solved by this invention is to provide a refined source apportionment method, system and storage medium for phosphorus pollution based on the coupled non-conservative transformation mechanism of phosphorus speciation and soluble organic matter, so as to solve the problem that the existing receptor model relies too much on the conservative assumption of tracers and lacks quantitative characterization means of phosphorus speciation mechanism, resulting in inaccurate source apportionment.

[0006] The technical solution adopted by this invention to solve its technical problem is: a refined source apportionment method for phosphorus pollution, comprising the following steps: Phosphorus speciation abundance data, driving factor data, and microbial characteristic data consisting of microbial community and phosphorus cycling functional genes were collected from samples from various receptor sampling points and potential pollution sources in the river channel, and phosphorus speciation abundance matrix, driving factor matrix, and microbial characteristic matrix were formed respectively; wherein, the driving factor data includes environmental factor data and relative abundance data of different degradation activity soluble organic matter components. The macroscopic source-sink effect coefficients of each phosphorus form were calculated using the aforementioned driving factor matrix; The biochemical effect coefficients of each phosphorus form were calculated using the aforementioned microbial characteristic matrix; The phosphorus speciation abundance matrix is ​​subjected to a central log-ratio transformation. Then, based on the transformation results, the nonconservative risk and environmental noise ratio of each phosphorus speciation are calculated. Effective tracers are selected according to a preset threshold, and the transformed abundance parameters of the effective tracers are extracted. The calculated macroscopic source-sink effect coefficients of each phosphorus form are coupled with the corresponding biochemical effect coefficients of the same phosphorus form to generate the comprehensive effect coefficients of each phosphorus form. For the selected effective tracers, the initial source spectrum mean and standard deviation of the corresponding pollution sources are synchronously scaled and corrected using the comprehensive effect coefficient, and then re-normalized to obtain the dynamic corrected mean and dynamic corrected standard deviation of each effective tracer, thereby reconstructing the dynamic corrected source spectrum matrix; and the equivalent fractionation coefficient of each pollution source at each sampling point and the corresponding hydrological season is extracted, wherein the equivalent fractionation coefficient is the difference between the dynamic corrected mean and the initial source spectrum mean. The abundance parameters, the dynamically corrected source spectrum matrix, and the equivalent fractionation coefficient are simultaneously input into the Bayesian mixture model. The iterative sampling of the posterior probability distribution is performed through its built-in iterative sampling algorithm. The contribution rate distribution of each pollution source to a specific phosphorus form in the receptor and the corresponding mean are output through mass-weighted calculation. Based on the mean of the posterior probability distribution of the contribution rate of each pollution source and the dynamically corrected source spectrum matrix, the theoretical predicted abundance of each phosphorus species in the receptor is calculated, and the basic residual is obtained by subtracting the theoretical predicted abundance from the actual observed value. A univariate linear regression model is constructed using the theoretically predicted abundance as the independent variable and the actual observed value as the dependent variable. The global analysis accuracy is quantified by the coefficient of determination and root mean square error of the model. The basic residuals are standardized and significance tested to identify significant outliers. If both the global analysis accuracy and the number of significant outliers meet their respective preset thresholds, the source apportionment result is deemed reliable, and the posterior probability distribution of the contribution rate of each pollution source is output. If either the global analysis accuracy or the number of significant outliers does not meet its respective preset threshold, the current source apportionment result is deemed unreliable, and the posterior probability distribution of the contribution rate of each pollution source is not output.

[0007] Furthermore, the phosphorus speciation in the phosphorus speciation abundance data includes two major categories: inorganic phosphorus and organic phosphorus. Inorganic phosphorus includes exchangeable phosphorus, aluminum-bound phosphorus, iron-bound phosphorus, and calcium-bound phosphorus; organic phosphorus includes sodium bicarbonate-extractable organic phosphorus, hydrochloric acid-extractable organic phosphorus, fulvic acid-bound organic phosphorus, humic acid-bound organic phosphorus, and residual organic phosphorus.

[0008] Furthermore, the step of calculating the macroscopic source-sink effect coefficients of each phosphorus form using the driving factor matrix includes: The feature importance score of each driving factor variable in the driving factor matrix to phosphorus form transformation is extracted using the lightweight gradient boosting machine algorithm. The correlation coefficient between each driving factor variable and phosphorus form transformation is calculated by combining the Mantell test. The feature importance score is multiplied by the correlation coefficient to obtain the comprehensive weight of each driving factor variable. The relative abundance of each degradation-active soluble organic matter component is multiplied by the combined weight of the driving factor variable corresponding to that component, and all products are summed to obtain the driving force score. The driving force score is normalized by range normalization to obtain the normalized driving force score. ; Calculate the macroeconomic source-sink effect coefficient based on the normalized driving force score. : ; In the formula: Source-sink direction discrimination coefficient; To correct the step size parameter.

[0009] Further, the step of calculating the biochemical effect coefficients of each phosphorus form using the microbial feature matrix includes: Calculate the global driving weights of microbial community and phosphorus cycling functional genes on each phosphorus form in the microbial feature matrix; The biochemical driving potential of the microbial community and phosphorus cycling functional genes was calculated, and the biochemical effect coefficient of each phosphorus form was obtained by combining the ecological efficiency score.

[0010] Further, the calculation of the global driving weights of microbial community and phosphorus cycling functional genes on each phosphorus form in the microbial feature matrix includes: Using the microbial community, phosphorus cycle functional genes, and phosphorus forms as latent variable nodes, the path relationship between the microbial community node and the functional gene node and the phosphorus form node is set. Using a partial least squares path model, the path coefficients from the microbial community nodes to the phosphorus morphology nodes and the path coefficients from the phosphorus cycle functional gene nodes to the phosphorus morphology nodes are extracted, and the sum of the absolute values ​​of the two types of path coefficients is used as the global driving weight.

[0011] Furthermore, the ecological efficiency score is obtained in the following way: A collinear network of phosphorus speciation, microbial community, and functional genes was constructed. The proportion of negatively correlated connections in the network was extracted, the structural stability efficiency was calculated, and then combined with the seasonal metabolic intensity index to obtain the ecological efficiency score.

[0012] Furthermore, the calculation of the biochemical driving potential of the microbial community and phosphorus cycling functional genes, combined with the eco-efficiency score, yields the biochemical effect coefficients of each phosphorus form, including: Using the random forest algorithm, with the relative abundance of microbial genera and the relative abundance of phosphorus cycling functional genes in the microbial community as input features and the phosphorus form transformation intensity as the prediction target, the feature contribution weights of each microbial genera and each phosphorus cycling functional gene to phosphorus form transformation are extracted. The correlation coefficient between the relative abundance of the phosphorus cycle functional genes and the phosphorus form conversion intensity was calculated by correlation analysis, and the driving potential of the functional genes was calculated based on the correlation coefficient. The microbial driving potential is calculated based on the correlation coefficient between the relative abundance at the microbial genus level and the phosphorus form transformation intensity, as well as the microbial component in the feature contribution weight. The microbial driving potential is multiplied by the functional gene driving potential and scaled to the 0.5-1.5 range using a linear mapping to obtain the standardized comprehensive biochemical driving potential. Multiply the global driving weight, the comprehensive biochemical driving potential, and the ecological efficiency score to obtain the biochemical incremental score. Then, multiply the biochemical incremental score by the biochemical scaling step size parameter and add 1 to obtain the biochemical effect coefficient of each phosphorus form.

[0013] Furthermore, the formula for calculating the functional gene drive potential is as follows: ; In the formula, The functional gene driving potential corresponding to the i-th phosphorus form, denoted as the correlation coefficient between the relative abundance of the j-th phosphorus cycling functional gene and the phosphorus form conversion intensity. Let be the average relative abundance of the j-th phosphorus cycle functional gene in all samples, and m be the total number of phosphorus cycle functional genes.

[0014] Furthermore, the formula for calculating the microbial-driven potential is as follows: ); In the formula, The microbial driving potential corresponding to the i-th phosphorus form, is the correlation coefficient between the relative abundance of the k-th microbial genus and the phosphorus form transformation intensity. Let l be the feature contribution weight of the k-th microbial genus to phosphorus form transformation obtained by the random forest algorithm, where l is the total number of microbial genera.

[0015] This invention also provides a refined source apportionment system for phosphorus pollution, used to implement the above-mentioned refined source apportionment method for phosphorus pollution. It is characterized by comprising: a data acquisition module, a macroscopic source-sink effect quantification module, a microscopic biochemical effect quantification module, a tracer screening module, a comprehensive effect coefficient generation module, a source spectrum dynamic correction module, a Bayesian contribution rate calculation module, and a result verification module.

[0016] The present invention also provides a computer-readable storage medium storing a computer program, which, when executed by a processor, implements the above-described method for refined source analysis of phosphorus pollution.

[0017] The beneficial effects of this invention are: Compared to the traditional source apportionment model's default assumption of tracer conservatism, this invention quantifies the intensity of non-conservative transformation of phosphorus speciation during migration by constructing a coupled model of macroscopic source-sink effects and microscopic biochemical effects. By introducing soluble organic matter molecular descriptors, microbial communities, and phosphorus cycle functional gene characteristics, it effectively solves the problem of apportionment bias caused by the inability to characterize non-conservative processes in complex biochemical environments.

[0018] This invention constructs a multidimensional environmental fingerprint database containing phosphorus speciation and dissolved organic matter molecules, significantly expanding the tracer information capacity compared to single isotope tracing. By utilizing a two-dimensional evaluation matrix based on the boundary crossing ratio and environmental noise ratio to perform tracer screening, redundant indicators severely affected by environmental fluctuations can be accurately eliminated. Experiments demonstrate that even in watershed scenarios with high land use heterogeneity and significant overlap of source spectral feature signals, this method maintains high spatial resolution and significantly reduces the identification ambiguity between pollution sources with similar contribution rates.

[0019] This invention utilizes the comprehensive effect coefficient to dynamically and mathematically reconstruct the static prior source spectrum of potential pollution sources, and quantitatively extracts the equivalent fractionation coefficients to input into the Bayesian mixture model. This solves the problem that traditional models fail to converge or output physically meaningless posterior probabilities when dealing with highly reactive phosphorus pollution scenarios due to the significant deviation of the measured signal at the receptor end from the range of the static source spectrum. Through mathematical reconstruction, the posterior output of the model not only has statistical significance, but also conforms to the actual physicochemical laws of phosphorus transport and transformation in the watershed. Attached Figure Description

[0020] Figure 1 This is a flowchart of the present invention; Figure 2 This is a regression analysis graph of measured and predicted values ​​in an embodiment of the present invention; Figure 3 This is a residual and outlier analysis chart from an embodiment of the present invention; Figure 4 This is a spatiotemporal distribution map of the contribution rate of farmland sources according to an embodiment of the present invention; Figure 5 This is a spatiotemporal distribution diagram of the contribution rate of human / engineering sources in an embodiment of the present invention; Figure 6 This is a spatiotemporal distribution diagram of the contribution rate of high organic matter sources in an embodiment of the present invention; Figure 7 This is a spatiotemporal distribution diagram of the sediment source contribution rate in an embodiment of the present invention. Detailed Implementation

[0021] The present invention will be further described below with reference to the accompanying drawings and embodiments.

[0022] like Figure 1 As shown, the refined source apportionment method for phosphorus pollution of the present invention includes the following steps: S1. Collect phosphorus speciation abundance data, driving factor data, and microbial characteristic data composed of microbial community and phosphorus cycling functional genes from samples of various receptor sampling points and potential pollution sources in the river channel, and form phosphorus speciation abundance matrix, driving factor matrix and microbial characteristic matrix respectively; wherein, the driving factor data includes environmental factor data and relative abundance data of different degradation activity soluble organic matter components. S2. Calculate the macroscopic source-sink effect coefficients of each phosphorus form using the driving factor matrix; S3. Calculate the biochemical effect coefficients of each phosphorus form using the aforementioned microbial characteristic matrix; S4. Perform a central log-ratio transformation on the phosphorus speciation abundance matrix, then calculate the nonconservative risk and environmental noise ratio of each phosphorus speciation based on the transformation result, select effective tracers according to a preset threshold, and extract the transformed abundance parameters of the effective tracers. S5. Couple the calculated macroscopic source-sink effect coefficients of each phosphorus form with the corresponding biochemical effect coefficients of the same phosphorus form to generate the comprehensive effect coefficients of each phosphorus form. S6. For the selected effective tracers, the initial source spectrum mean and standard deviation of the corresponding pollution sources are synchronously scaled and corrected using the comprehensive effect coefficient, and then re-normalized to obtain the dynamic corrected mean and dynamic corrected standard deviation of each effective tracer, thereby reconstructing the dynamic corrected source spectrum matrix; and extracting the equivalent fractionation coefficient of each pollution source at each sampling point and the corresponding hydrological season, wherein the equivalent fractionation coefficient is the difference between the dynamic corrected mean and the initial source spectrum mean; S7. The abundance parameters, the dynamically corrected source spectrum matrix and the equivalent fractionation coefficient are simultaneously input into the Bayesian mixture model. The iterative sampling of the posterior probability distribution is performed through its built-in iterative sampling algorithm. The contribution rate distribution of each pollution source to the specific phosphorus form in the receptor and the corresponding mean are output through mass weighting calculation. S8. Based on the mean of the posterior probability distribution of the contribution rate of each pollution source and the dynamically corrected source spectrum matrix, calculate the theoretical predicted abundance of each phosphorus species in the receptor, and obtain the basic residual by subtracting the theoretical predicted abundance from the actual observed value. S9. Construct a univariate linear regression model using the theoretically predicted abundance as the independent variable and the actual observed value as the dependent variable. Quantify the global analysis accuracy using the model's coefficient of determination and root mean square error. Standardize and perform significance testing on the basic residuals to identify significant outliers. If both the global analysis accuracy and the number of significant outliers meet their respective preset thresholds, the source apportionment result is deemed reliable, and the posterior probability distribution of the contribution rate of each pollution source is output, along with a conclusion that the verification has passed. If either the global analysis accuracy or the number of significant outliers does not meet its respective preset threshold, the current source apportionment result is deemed unreliable.

[0023] If the current source resolution result is determined to be unreliable, the following troubleshooting mechanism can be executed: Low global resolution accuracy or excessive local anomalies are usually due to the introduction of sudden, sporadic point sources or specific micro-pollutants not covered by the model into the recipient water body, or extreme environmental physical noise and strong non-conservative behavior at individual sampling points, causing the actual observed abundance to deviate significantly from the dynamic source spectrum optimization range reconstructed by the multi-mechanism coupling model. In cases where verification fails, the specific handling steps are as follows: (1) Automatically suspend the output of the current source analysis result and prompt the administrator that there is a risk of analysis uncertainty at the current location; (2) Collect water samples from the corresponding abnormal sampling points or key nodes along the route during this period and seal them physically for subsequent higher-precision laboratory physicochemical and isotope re-analysis. (3) Based on the abnormal locations indicated by the alarm and the subsequent high-precision analysis results, the newly identified specific pollution source feature fingerprints are incorporated into the potential pollution source database to dynamically correct the prior source spectrum boundary.

[0024] This invention collects phosphorus speciation abundance data, driving factor data, and microbial characteristic data and forms matrices for each. The driving factor matrix is ​​used to quantify macroscopic source-sink effect coefficients, and the microbial characteristic matrix is ​​used to quantify biochemical effect coefficients. These two are coupled to generate a comprehensive effect coefficient, which dynamically corrects the initial source spectrum. Simultaneously, effective tracers are screened and equivalent fractionation coefficients are extracted through central logarithmic ratio transformation, non-conservative risk, and environmental noise ratio. A Bayesian mixture model and iterative sampling are combined to output the pollution source contribution rate distribution. Finally, a regression model and residual significance test are used to achieve global accuracy quantification and local outlier identification. This forms a complete technical chain from data acquisition, mechanism coupling, dynamic correction to multi-dimensional verification, effectively overcoming the shortcomings of traditional methods that rely on conservative assumptions and lack quantitative characterization of transformation mechanisms, significantly improving the accuracy and reliability of phosphorus pollution source apportionment.

[0025] The following provides a detailed explanation of each of the above steps.

[0026] Detailed explanation of step S1: The phosphorus speciation data can be categorized into two main types: inorganic phosphorus and organic phosphorus. Inorganic phosphorus includes exchangeable phosphorus (Ex-P), aluminum-bound phosphorus (Al-P), iron-bound phosphorus (Fe-P), and calcium-bound phosphorus (Ca-P); and organic phosphorus extracted with sodium bicarbonate (NaHCO3-P), organic phosphorus extracted with hydrochloric acid (HCl-P), organic phosphorus bound with fulvic acid (Ful-P), organic phosphorus bound with humic acid (Hum-P), and organic phosphorus in residue form (Res-P). Among these, exchangeable phosphorus (Ex-P), aluminum-bound phosphorus (Al-P), and iron-bound phosphorus (Fe-P) are active phosphorus speciations, while the rest are inactive phosphorus speciations.

[0027] In this embodiment of the invention, the relative abundance of exchangeable phosphorus, aluminum-bound phosphorus, iron-bound phosphorus, and calcium-bound phosphorus in the sample is obtained by graded continuous extraction method, and the relative abundance of sodium bicarbonate-extractable organic phosphorus, hydrochloric acid-extractable organic phosphorus, fulvic acid-bound organic phosphorus, humic acid-bound organic phosphorus, and residual phosphorus is obtained by calculating the difference between total phosphorus and inorganic phosphorus.

[0028] In this invention, environmental factors refer to physical and chemical parameters that characterize the water environment of the recipient water body and the pollution source end. These parameters, together with dissolved organic matter components, constitute driving factors, including but not limited to water temperature, pH value, dissolved oxygen, conductivity, turbidity, flow velocity, flow rate, permanganate index, ammonia nitrogen, total nitrogen, total phosphorus, redox potential, and chlorophyll a. The environmental factors listed above cover the conventional core indicators in water environment quality assessment. In actual implementation, other factors (such as nitrate nitrogen, phosphate, heavy metal ions, etc.) can be selected according to the characteristics of watershed pollution and monitoring conditions. As long as they can participate in the subsequent comprehensive weighting and macro-source-sink effect coefficient calculation, they should be considered as environmental factors covered by this invention.

[0029] The relative abundance data of the different degradation-active soluble organic matter components can be obtained by pure water extraction combined with Fourier transform ion cyclotron resonance mass spectrometry analysis. Specifically, the relative abundance of nominal carbon oxidation state, aromaticity index and characteristic compound components in each sample is calculated by high-resolution mass spectra, thereby constructing a molecular characteristic parameter matrix of soluble organic matter.

[0030] The microbial feature matrix was obtained through high-throughput sequencing and metagenomic sequencing. High-throughput sequencing was used to obtain the distribution structure of the genus-level microbial community; metagenomic sequencing was used to extract the relative abundance of functional genes related to inorganic phosphorus dissolution and organic phosphorus mineralization, thus forming a microbial feature matrix composed of microbial community and phosphorus cycling functional genes.

[0031] In the phosphorus speciation abundance matrix, each row corresponds to a sample (including samples from pollution sources and receptors), and each column corresponds to a phosphorus speciation (such as exchangeable phosphorus, aluminum-bound phosphorus, iron-bound phosphorus, calcium-bound phosphorus, sodium bicarbonate-extractable organic phosphorus, hydrochloric acid-extractable organic phosphorus, fulvic acid-bound organic phosphorus, humic acid-bound organic phosphorus, and residual organic phosphorus). The matrix element values ​​are the relative abundance of the corresponding phosphorus speciation in that sample. In the driving factor matrix, each row corresponds to a sample (in the same row order as the phosphorus speciation abundance matrix), and each column corresponds to a driving factor variable (including water temperature, pH, dissolved oxygen, conductivity, turbidity, flow rate, permanganate index, ammonia nitrogen, total nitrogen, total phosphorus, and redox potential). The matrix contains environmental factors such as chlorophyll a, and the relative abundance of soluble organic matter components with different degradation activities, such as aliphatic compounds, high-oxygen phenols, low-oxygen phenols, and polycyclic aromatic hydrocarbons. The matrix element values ​​are the measured values ​​of the corresponding environmental factors or the relative abundance of soluble organic matter components in the sample. In the microbial feature matrix, each row corresponds to one sample (in the same order as the rows in the matrix above), and each column corresponds to one microbial feature variable (including the relative abundance at the genus level, and the relative abundance of phosphorus cycle functional genes such as inorganic phosphorus dissolving genes, organic phosphorus mineralization genes, phosphorus transport genes, phosphorus regulation genes, and carbohydrate active enzyme genes). The matrix element values ​​are the relative abundance of the corresponding genus or functional gene in the sample.

[0032] Detailed explanation of step S2: In this embodiment of the invention, step S2, the step of calculating the macroscopic source-sink effect coefficients of each phosphorus form using the driving factor matrix, specifically includes: S21. Using the lightweight gradient boosting machine algorithm, extract the feature importance scores of each driving factor variable in the driving factor matrix for phosphorus speciation conversion. Combine this with the Mantell test to calculate the correlation coefficient between each driving factor variable and phosphorus speciation conversion. Multiply the feature importance scores by the correlation coefficients to obtain the comprehensive weights of each driving factor variable. ; In the formula, The combined driving weight of the j-th driving factor variable on the i-th phosphorus form; The feature importance score of the j-th driving factor variable to the i-th phosphorus form, obtained by quantification through the LGBM model; The coupling index is the correlation degree r value of a significant correlation (p < 0.05) in the Mantel test. If it is not significant, it is taken as 0. For driving factors with p values ​​close to the significance level (e.g., 0.05 < p < 0.1), their contribution can be retained as appropriate.

[0033] S22. Multiply the relative abundance of each degradation active soluble organic matter component by the combined weight of the driving factor variable corresponding to that component, and sum all products to obtain the driving force score: ; In the formula, The original driving force score for the i-th phosphorus form in the k-th sample; denoted as , where is the relative abundance of the j-th soluble organic matter component in the k-th sample; n is the total number of driving factor variables.

[0034] S23. Perform range normalization on the driving force score to obtain the normalized driving force score. : ; In the formula, and Let be the minimum and maximum values ​​of the driving force score for the i-th phosphorus form among all samples in a specific season. = hour, A value of 0.5 indicates that the driving force of this phosphorus form has no significant spatial variation, and its macroscopic source-sink effect is determined only by the source-sink direction discrimination coefficient. Then, the macroeconomic source-sink effect coefficient is calculated based on the normalized driving force score. : ; In the formula: The source-sink direction discrimination coefficient is determined by the positive or negative value of the relative change rate of the i-th phosphorus form in the k-th sample: 1 when it is greater than 0, -1 when it is less than 0, and 0 when it is equal to 0 or missing. To correct the step size parameter, the value range is 0.1 to 1.0, and the preferred value is 0.5.

[0035] To prevent outliers from causing distortion in the source spectrum correction, i and k are truncated, that is The value range is constrained to [0.1, 2.0].

[0036] In some embodiments, the above normalization process can also be replaced by Z-score standardization (i.e., subtracting the mean and dividing by the standard deviation) or decimal scaling normalization instead of range normalization; the calculation of the driving force score can also be achieved by directly using weighted summation without explicit normalization, and a similar effect can be achieved by adjusting the correction step size η. In addition, the calculation of the comprehensive weight can also use random forest feature importance or Shapley value instead of the LightGBM algorithm, as long as the contribution intensity of each driving factor to phosphorus speciation can be quantified.

[0037] Detailed explanation of step S3: In this embodiment of the invention, step S3: calculating the biochemical effect coefficients of each phosphorus form using the microbial feature matrix specifically includes: S31. Calculate the global driving weights of microbial community and phosphorus cycling functional genes on each phosphorus form in the microbial feature matrix. S32. Calculate the biochemical driving potential of the microbial community and phosphorus cycling functional genes, and combine the ecological efficiency score to obtain the biochemical effect coefficient of each phosphorus form.

[0038] To more scientifically integrate the comprehensive structural information and functional gene abundance information of the microbial community, and to provide robust and interpretable weight parameters for subsequent quantification of microscopic biochemical effects, this embodiment of the invention uses the microbial community, phosphorus cycling functional genes, and phosphorus speciation as latent variable nodes, and sets the path relationships between microbial community nodes and functional gene nodes pointing to phosphorus speciation nodes: the latent variables of the microbial community and phosphorus cycling functional genes are set as exogenous variables, and the latent variables of phosphorus speciation (divided into active phosphorus and inactive phosphorus) are set as endogenous variables. Using a partial least squares path model (PLS-PM), the path coefficients from the microbial community nodes to the phosphorus speciation nodes and from the phosphorus cycling functional gene nodes to the phosphorus speciation nodes are extracted, and the sum of the absolute values ​​of the two types of path coefficients is used as the global driving weight. ; In the formula, The global weight of the microbial community on the c-th phosphorus form (active or inactive phosphorus); The path coefficients from microbial community nodes to the c-th type of phosphorus speciation node in the PLS-PM model; The significant path coefficients from functional gene nodes to the c-th type of phosphorus morphology node are given.

[0039] In this embodiment of the invention, step S32 specifically includes: S321. Using the random forest algorithm, with the relative abundance of microbial genera and the relative abundance of phosphorus cycling functional genes in the microbial feature matrix as input features, and with the phosphorus form transformation intensity as the prediction target, extract the feature contribution weights of each microbial genera and each phosphorus cycling functional gene to each phosphorus form transformation. S322. Calculate the correlation coefficient between the relative abundance of the phosphorus cycling functional genes and the phosphorus form conversion intensity through correlation analysis, and calculate the functional gene driving potential based on the correlation coefficient; specifically, the formula for calculating the functional gene driving potential is as follows: ; In the formula, The functional gene driving potential corresponding to the i-th phosphorus form, denoted as the correlation coefficient between the relative abundance of the j-th phosphorus cycling functional gene and the phosphorus form conversion intensity. Let be the average relative abundance of the j-th phosphorus cycle functional gene in all samples, and m be the total number of phosphorus cycle functional genes. S323. Calculate the microbial driving potential based on the correlation coefficient between the relative abundance at the microbial genus level and the phosphorus form transformation intensity, and the microbial component in the feature contribution weight; the formula for calculating the microbial driving potential is as follows: ); In the formula, The microbial driving potential corresponding to the i-th phosphorus form, is the correlation coefficient between the relative abundance of the k-th microbial genus and the phosphorus form transformation intensity. Let l be the feature contribution weight of the k-th microbial genus to phosphorus form transformation obtained by the random forest algorithm, where l is the total number of microbial genera. S324. Multiply the microbial driving potential by the functional gene driving potential, and scale to the range of 0.5–1.5 using a linear mapping to obtain the standardized comprehensive biochemical driving potential; if If the linear mapping result exceeds the 0.5~1.5 range, the boundary value of the range is directly taken to eliminate extreme anomalies in biochemical driving potential. ; In the formula, This represents the potential for a standardized, integrated biochemical process. It is a linear mapping function.

[0040] S325, Set the global driving weights The aforementioned integrated biochemical driving potential With the aforementioned ecological efficiency score Multiply to obtain the biochemical increment score. Then, the biochemical increment score is calculated. Multiply by the biochemical scaling step parameter Adding 1 gives the biochemical effect coefficients for each phosphorus form: ; ; In the formula, Let be the biochemical increment score of the i-th phosphorus form in season s, where The value of depends on whether the i-th phosphorus form belongs to the category of active or inactive phosphorus. The final output is the biochemical effect coefficient, which will be compared with the macroscopic source-sink effect coefficient. The common input is added to the multi-mechanism coupling model; γ is the biochemical scaling step size parameter, which is used to control the correction intensity of the biochemical effect coefficient to the source spectrum. The preferred value is 1.0, and the value range is 0.5 to 1.5.

[0041] In some embodiments, methods such as XGBoost, LightGBM, Gradient Boosting Tree (GBDT), LASSO regression, Elastic Net, Mutual Information, Distance Correlation Coefficient (dCor), Permutation Importance, Boruta algorithm, or Shapley Additive Explanations are used to replace the random forest algorithm in order to extract the characteristic contribution weights of each microbial genus and each phosphorus cycling functional gene to phosphorus form transformation.

[0042] In this embodiment of the invention, the ecological efficiency score is obtained in the following way: A collinear network of phosphorus speciation, microbial community, and functional genes was constructed. The proportion of negatively correlated connections in the network was extracted, and the structural stability efficiency was calculated. This efficiency was then combined with the seasonal metabolic intensity index to obtain the ecological efficiency score. The specific formula for the ecological efficiency score is as follows: ; ; In the formula, For the structural stability efficiency of the microbial network in season s, when When ≤0, take 0.01 or Take percentages, with values ​​ranging from 0 to 100; The proportion of negatively correlated connections in the co-occurrence network during season s; The final ecological efficiency score is derived by comprehensively considering seasonal characteristics; It is a seasonal metabolic intensity index used to reflect the promoting or inhibiting effects of environmental factors such as water temperature and flow rate on biological metabolism.

[0043] The following are The calculation process for the value is illustrated with an example.

[0044] First, the effect of temperature on microbial activity was calculated using the microbial metabolic temperature correction formula: ; Where θ is the temperature correction coefficient, taken as an empirical constant of 1.05; T is the average water temperature in a specific season; The average annual water temperature of the basin.

[0045] Taking a typical small watershed in the Sichuan Basin as an example, the average annual water temperature is about 18 degrees Celsius, and the average water temperature during the rainy season is usually between 23 and 25 degrees Celsius. Taking the median value of 24 degrees Celsius for calculation, The average water temperature during the dry season is usually between 11 and 13 degrees Celsius; taking the median of 12 degrees Celsius for calculation, .

[0046] Furthermore, considering the lower river flow velocity during the dry season, the extended interaction time between microorganisms and phosphorus forms can partially compensate for the inhibition of metabolic activity caused by low temperatures. To quantitatively determine the degree of compensation, the hydraulic residence time ratio method can be used: according to watershed hydrological monitoring data, the average flow velocity during the dry season is about 80% lower than that during the rainy season, and the hydraulic residence time is extended to five times that of the rainy season. Considering the diminishing marginal returns of the extended time effect, the compensation coefficient is set to the power of 0.4 of the residence time ratio (i.e., 5). 0.4 (≈1.9), and to avoid overcompensation, the upper limit of the compensation coefficient is set at 1.2. Therefore, the compensation coefficient after dry season compensation is calculated. : .

[0047] During the rainy season, the flow velocity is high and the temperature is suitable, so no additional compensation is needed. Take the calculated value of 1.27 directly.

[0048] Detailed explanation of step S4: Since the sum of the relative abundance data of phosphorus species is a constant (i.e., closure effect), directly using the raw abundance data will lead to multicollinearity. Therefore, a central log-ratio (CLR) transformation is performed on the raw relative abundance of phosphorus species in each sample to eliminate data closure constraints. The transformation formula is as follows: ; In the formula, For the first i The values ​​of phosphorus species after central logarithmic ratio transformation; denoted as the original relative abundance of this phosphorus form; D represents the total number of phosphorus forms.

[0049] Non-conservative risk quantification: Non-conservative risk measures the degree to which the distribution of a certain phosphorus form at the receptor end deviates from the coverage range of the pollution source fingerprint. The tracer effectiveness of this indicator is assessed by calculating the proportion of receptor samples falling outside the effective range (i.e., the out-of-bounds proportion). Specifically, the minimum and maximum values ​​of the specific phosphorus form transformation values ​​in each potential pollution source are first extracted, and an effective tracer range is constructed based on this range plus a 10% tolerance threshold. If a receptor sample falls outside this range, it is considered a significant transformation exceeding the model's correction capability. The out-of-bounds proportion calculation formula is as follows: ; In the formula, For the firsti The risk of non-conservatism in phosphorus species (i.e., the proportion of species crossing the boundary). The smaller the value, the lower the degree of conversion of the phosphorus form from the source to the acceptor, and the more suitable it is as a tracer; The number of samples in the receptor sample whose transformation value of the i-th phosphorus form deviates from the effective range is the source extreme value ± 10% range. The 10% is the default tolerance parameter, which can be adjusted within the range of 5% to 20% according to the hydrological variability of the watershed. This represents the total number of receptor samples.

[0050] Environmental noise ratio quantification: The environmental noise ratio is used to assess the relative magnitude of receptor-side variability to the inherent variability at the source end, preventing the source signal from being masked by environmental fluctuations. The calculation formula is: ; In the formula, For the first i The environmental noise ratio of phosphorus species; and The receptor sample set and the pollution source sample set are respectively the first two sets. i Standard deviation of phosphorus speciation values. The closer the value is to 1, the more comparable the receptor variation is to the source end; if it is too large (exceeding the preset threshold of 1.2), the receptor variation is much greater than the source end difference, and this indicator is not suitable as a tracer.

[0051] Two-dimensional filtering and parameter extraction: combined and Two core quantitative indicators are used to construct a two-dimensional evaluation matrix. Phosphorus speciation that simultaneously meets both indicators and is controlled within preset thresholds is identified as an effective tracer and included in subsequent source apportionment calculations. After screening, the CLR transform values ​​of effective tracers are extracted to obtain the transformed abundance parameters, which are used for subsequent source spectrum dynamic correction and Bayesian mixture model input.

[0052] In a preferred embodiment, the threshold for non-conservative risk can be set to 20%, based on the following: Bayesian mixture models are based on the law of mass conservation, and the receptor signal must fall within the convex hull formed by the pollution sources to be decomposed into a non-negative linear combination of the contribution rates of each source; if more than 20% of the receptor samples deviate from the convex hull boundary, the tracer conversion intensity exceeds the inherent characteristics of the source spectrum, causing the MCMC chain to get stuck at the boundary and the model to fail to converge. The environmental noise ratio threshold can be set to 1.2, based on the following: According to the error propagation law, when the environmental noise ratio threshold is 1.2, the extra variance at the receptor end has reached 44% of the inherent variance of the source spectrum, and the signal-to-noise ratio is too low, so the characteristic signal will be submerged. In practical applications, the above default values ​​can be adjusted according to the watershed heterogeneity and analytical accuracy requirements.

[0053] Detailed explanation of step S6: For the selected effective tracers, the combined effect coefficient is used. The initial source spectrum mean and standard deviation of the corresponding pollution sources are synchronously scaled and corrected to obtain the effective mean and effective standard deviation. The comprehensive effect coefficient... Generated in step S5, its subscript i represents the phosphorus speciation type, k represents the sampling scenario (sampling location and spatiotemporal combination), and s represents the hydrological season (rainy season or dry season), comprehensively reflecting the combined influence of macro-environmental driving forces and micro-biological driving forces on phosphorus speciation migration and transformation. The specific correction formula is as follows: ; ; In the formula, and These are the effective mean and effective standard deviation of the i-th phosphorus form in the v-th pollution source under the k-th sample scenario after migration and transformation; and These represent the initial mean and initial standard deviation of the i-th phosphorus form in the v-th pollution source, as determined by preliminary field surveys. The comprehensive effect coefficient of the i-th phosphorus form in the k-th sampling scenario and the s-th hydrological season is obtained by coupling the macro-source-sink effect coefficient and the biochemical effect coefficient (i.e., ).

[0054] Due to the above effective mean The mass balance constraint (the sum of the proportions of each phosphorus form should be 1) is not yet met, and renormalization is required. Simultaneously, to maintain the variability of the source spectrum, the effective standard deviation is equivalently transformed by the same proportion: ; ; In the formula, The mean of the i-th phosphorus form is dynamically corrected and finally input into the Bayesian source analysis model after renormalization. denoted as the dynamically corrected standard deviation for the i-th phosphorus species after equivalent normalization; D represents the total number of phosphorus species. The dynamically corrected mean and standard deviation of all effective tracers together constitute the dynamically corrected source spectrum matrix.

[0055] To adapt to the input structure of the Bayesian mixture model, the source spectrum shift caused by the non-conservative transformation is quantized into equivalent fractionation coefficients. The difference between the dynamically corrected mean and the initial source spectrum mean is calculated: ; In the formula, Let be the average equivalent fractionation coefficient of the i-th phosphorus form in the v-th pollution source under the k-th scenario. Finally... and and Together, these constitute the core prior parameters, which include macroscopic source-sink responses and microscopic biochemical effects, and are input into the subsequently constructed Bayesian mixture model.

[0056] Detailed explanation of step S7: In this embodiment of the invention, the Bayesian mixture model adopts the MixSIAR model framework and extends it to accept the equivalent fractionation coefficient as the source spectrum correction parameter.

[0057] Model Setup and Convergence Diagnosis: The model is run independently for each sampling point (single-point mode), using the Markov Chain Monte Carlo (MCMC) algorithm to iteratively sample the posterior probability distribution of each pollution source. The Gelman-Rubin diagnostic statistic (R) is extracted to evaluate the convergence of the MCMC chain, requiring that the R values ​​of all simulation parameters be strictly less than 1.05. If the diagnostic value exceeds this threshold, the model is considered not fully converged, and the iteration chain length is automatically increased until the convergence requirement is fully met.

[0058] Posterior probability distribution calculation: The final contribution rate of a specific phosphorus species in the acceptor is affected by the overall mixing ratio of the pollution sources, and also strictly depends on the abundance of that species in the corresponding pollution source. Therefore, after extracting the source mixing ratio generated in each simulation in the MCMC chain, it is necessary to perform mass weighting and re-normalization calculations in conjunction with the effective abundance of the source spectrum. The calculation formula is: ; ; In the formula, The overall mixing ratio of the v-th pollution source is calculated in the m-th iteration of the MCMC chain. is the effective abundance parameter of the i-th phosphorus form in the v-th pollution source; The weight of the absolute mass contribution of the v-th pollution source to the i-th phosphorus form in the m-th iteration; is the final contribution rate of the v-th pollution source to the i-th phosphorus form in a specific receptor after normalization; V is the total number of pollution sources input into the model.

[0059] Posterior distribution statistics and output: The set of contribution rates calculated based on the above iterative calculations. The posterior probability distribution matrix of the contribution rate of each pollution source to a specific phosphorus form at each sampling point is obtained. Through statistical calculation, the mean, standard deviation, and a full set of characteristic quantiles (including minimum, 2.5%, 5%, 25%, median, 75%, 95%, 97.5%, and maximum) of this probability distribution are extracted, thereby accurately determining the absolute level of the pollution source's contribution and the uncertainty boundary of the model prediction.

[0060] Detailed explanation of step S8: Based on the mean of the posterior probability distribution of the contribution rates of each pollution source and the dynamically corrected source spectrum matrix, the theoretical predicted abundance of each phosphorus species in the receptor is calculated. The specific formula is: ; In the formula, The theoretically predicted relative abundance of the i-th phosphorus species in the k-th sample; The total contribution rate of the vth pollution source to the kth sample, calculated by the model; is the effective abundance parameter of the i-th phosphorus form in the v-th pollution source; This represents the total number of pollution sources by category.

[0061] Then, the difference between this and the observed values ​​obtained from the actual test is calculated to obtain the basic residual. ; In the formula, The basic residual of the i-th phosphorus form in the k-th sample; This represents the relative abundance of receptor phosphorus species as observed in practice.

[0062] Detailed explanation of step S9: Global resolution precision quantification: based on the theoretically predicted relative abundance of receptor phosphorus species The relative abundance of receptor phosphorus species is the actual observed value. Using the rainy season and dry season datasets as the dependent variable, construct univariate linear regression models separately. Extract the coefficient of determination. The two metrics, along with the root mean square error (RMSE), are used to evaluate the model's ability to explain the global data and the spatial distribution variation of individual phosphorus speciation.

[0063] In a preferred embodiment, the determination coefficient The threshold can be set to 0.85, based on the fact that highly reactive non-conservative elements are strongly affected by biogeochemical cycles, and the industry-recognized convergence standard is... ≥0.80; This invention is based on dynamic correction of the mechanism, raising the standard to 0.85 to demonstrate that the dynamic correction mechanism significantly improves the fitting accuracy.

[0064] Standardized residual score calculation and local outlier identification: To accurately identify sporadic transformation events or specific, minute pollution sources within the recipient environment, the basic residuals are standardized based on specific hydrological seasons and phosphorus speciation, and their statistical significance level deviating from model predictions is calculated. The calculation formula is as follows: ; ); In the formula: The standardized residual score; and The following are the subsections for each hydrological season (rainy season or dry season). i Mean and standard deviation of the basic residuals of phosphorus speciation; This represents the two-sided significance level test value corresponding to the residual score. F This is the cumulative distribution function of the standard normal distribution.

[0065] Based on the above calculation results, a residual distribution matrix is ​​constructed, and the basic residuals are... The negative logarithm of the significance level is used as the horizontal axis. Used as the ordinate. Set P < 0.01 (i.e., -). >2) is the significance threshold. Abnormal sampling points exceeding this threshold are extracted, thereby identifying receptor samples in the spatial distribution that may have analytical biases. Combining the global analytical accuracy and the results of local anomaly identification, if both the global analytical accuracy and the number of significant anomalies meet their respective preset thresholds, the source apportionment result is deemed reliable, and the posterior probability distribution of the contribution rate of each pollution source is output, along with a conclusion that the verification has been passed. If either the global analytical accuracy or the number of significant anomalies does not meet its respective preset threshold, the current source apportionment result is deemed unreliable, and the posterior probability distribution of the contribution rate of each pollution source is not output.

[0066] In a preferred embodiment, the threshold for the proportion of significant outliers can be set to 5%, based on the following: when P<0.01 is set as the criterion for determining significant outliers, even if the model is completely correct, about 1% of the samples will still be misclassified as outliers due to random fluctuations; relaxing the tolerance to 5% corresponds to the random misclassification rate at the statistical significance level α=0.05, which can effectively identify real outliers and avoid oversensitivity.

[0067] This invention also provides a refined source apportionment system for phosphorus pollution, used to implement the above-mentioned refined source apportionment method for phosphorus pollution, comprising: The data acquisition module is configured to perform step S1; The macro-source-sink effect quantification module is configured to execute step S2; A module for quantifying microscopic biochemical effects is configured to perform step S3; A tracer screening module is configured to perform step S4; A comprehensive effect coefficient generation module is configured to execute step S5; The source spectrum dynamic correction module is configured to execute step S6; The Bayesian contribution rate calculation module is configured to perform step S7; and, The result verification module is configured to execute steps S8 and S9.

[0068] The present invention also provides a computer-readable storage medium storing a computer program, which, when executed by a processor, implements the above-described method for refined source analysis of phosphorus pollution.

[0069] Example: Taking a typical agricultural watershed in the Sichuan Basin as an example, the method described in this invention was used to perform refined source apportionment of phosphorus pollution.

[0070] Unspoiled samples were collected from river receptor sampling points and surrounding potential pollution sources within the watershed during both the rainy and dry seasons. Potential pollution sources were categorized into four groups: high organic matter group, agricultural source group, human activity input group, and endogenous sediment group.

[0071] The relative abundance of inorganic phosphorus forms such as exchangeable phosphorus (Ex-P), aluminum-bound phosphorus (Al-P), iron-bound phosphorus (Fe-P), and calcium-bound phosphorus (Ca-P) in the samples was determined by a fractional continuous extraction method, as well as organic phosphorus forms such as sodium bicarbonate-extractable organic phosphorus (NaHCO3-P), hydrochloric acid-extractable organic phosphorus (HCl-P), fulvic acid-bound organic phosphorus (Ful-P), humic acid-bound organic phosphorus (Hum-P), and residual organic phosphorus (Res-P) in the samples, forming a phosphorus form abundance matrix.

[0072] Fourier transform ion cyclotron resonance mass spectrometry (FT-ICR-MS) was used to obtain the nominal carbon oxidation state, aromaticity index, and relative abundance of different degradation active components (aliphatic compounds, high-oxygen / low-oxygen phenols, polycyclic aromatic hydrocarbons, etc.) of dissolved organic matter (DOM), forming a DOM molecular characteristic parameter matrix (i.e., the organic matter component part of the driving factor matrix). Simultaneously, environmental factors (water temperature, pH, dissolved oxygen, conductivity, turbidity, flow rate, permanganate index, ammonia nitrogen, total nitrogen, total phosphorus, etc.) data were collected, collectively forming the driving factor matrix.

[0073] High-throughput sequencing and metagenomic sequencing were used to obtain the relative abundance of genus-level microbial community structure and phosphorus cycle functional genes (inorganic phosphorus dissolution, organic phosphorus mineralization, phosphorus transport, etc.) and construct a microbial characteristic matrix.

[0074] The Lightweight Gradient Boosting Machine (LightGBM) algorithm was used to extract the characteristic importance scores of each driving factor (environmental factors and DOM components) on phosphorus speciation, and the comprehensive weight was calculated using the Mantel test. Driving force scores were calculated based on the relative abundance of different degradation-active DOM components and the comprehensive weight. After range normalization and combining with the source-sink direction discrimination coefficient (valued at ±1 or 0 based on the relative change rate of phosphorus speciation), the macroscopic source-sink effect coefficients of each phosphorus speciation were obtained. (Correction step size η=0.5, The value is constrained to be between 0.1 and 2.0.

[0075] Using microbial community, phosphorus cycle functional genes, and phosphorus form as latent variable nodes, a partial least squares path model (PLS-PM) was used to extract path coefficients and calculate global driving weights. (Active phosphorus form is 0.64).

[0076] The random forest algorithm was used to extract the feature contribution weights at the genus level and for functional genes, and the microbial driving potential and functional gene driving potential were calculated separately. After multiplication, the normalized comprehensive biochemical driving potential was obtained by linear mapping. (Scale to 0.5~1.5).

[0077] Construct a collinear network and calculate the structural stability efficiency based on the proportion of negatively correlated connections. Combined with seasonal metabolic intensity index (1.27 in the rainy season, 0.85 in the dry season) to obtain the ecological efficiency score (0.549 in the rainy season, 0.265 in the dry season). The biochemical effect coefficient is finally calculated using the following formula. Where γ = 1.0.

[0078] The phosphorus speciation abundance matrix was subjected to a central log-ratio (CLR) transformation to calculate the risk of nonconservatism. (Boundary crossing ratio) and environmental noise ratio .by 20% and Using a threshold value, five effective tracers—Ex-P, Al-P, Fe-P, Ca-P, and Res-P—were selected during the rainy season, while four effective tracers—Fe-P, Res-P, NaHCO3-P, and Hum-P—were selected during the dry season. The abundance parameters of the effective tracers after CLR transformation were extracted.

[0079] The comprehensive effect coefficient is obtained by coupling the macroeconomic source-sink effect coefficient with the biochemical effect coefficient. In this embodiment, the comprehensive effect coefficient of Ex-P during the rainy season is 1.843.

[0080] The initial source spectrum mean and standard deviation are synchronously scaled and corrected, then renormalized by 1.27 to obtain the dynamically corrected source spectrum matrix, and the equivalent fractionation coefficients are extracted. In this embodiment, the equivalent fractionation coefficient of Ex-P in the farmland source group during the rainy season is 0.042.

[0081] The transformed abundance parameters, dynamically corrected source spectrum matrix, and equivalent fractionation coefficients of the effective tracer were input into the MixSIAR model. A single-point model was used to independently analyze each sampling point, and the MCMC algorithm (3 chains, 100,000 iterations per chain, 50,000 preheating iterations) was used for posterior probability sampling. The Gelman-Rubin statistic remained below 1.03, strictly meeting the convergence threshold of less than 1.05, indicating model convergence. The posterior distribution (mean and 95% confidence interval) of the contribution rate of each pollution source to different phosphorus forms was output through mass-weighted calculation.

[0082] The theoretically predicted abundance is calculated based on the posterior contribution rate mean and dynamically corrected source spectrum parameters, and the difference between this and the actual observed values ​​is used to obtain the basic residual. A univariate linear regression model is constructed (predicted values ​​are independent variables, observed values ​​are dependent variables) for rainy and dry seasons. All values ​​were greater than 0.85, and the RMSE was less than 0.12, indicating that the model has high global analytical accuracy. After standardizing the basic residuals, P < 0.01 (i.e., - >2) Local outliers are identified using thresholding. The results of theoretical abundance regression and local outlier analysis in this embodiment are shown below. Figure 2 , Figure 3 No significant anomalies were found in the figure, verifying the reliability of the source resolution results.

[0083] Cross-validation: The pollution source contribution rate of each sampling point output by the model was extracted, and a spatial distribution heatmap of the contribution rate of each type of pollution source was generated using a spatial interpolation algorithm. This spatial distribution heatmap was then spatially topologically overlaid with high-precision land use classification data and the coordinates of the pollution source field survey. The cross-validation results of the spatial distribution of pollution source contribution rates are shown below. Figures 4 to 7 The results show that the core distribution areas of high contribution rates of each pollution source have a high degree of spatial consistency with the corresponding land use types (such as the distribution of high-value areas of farmland sources and cultivated land, and the distribution of high-value areas of high-organic-matter sources and rural settlements), thus cross-validating the reliability of the source apportionment results of the method of the present invention.

Claims

1. A refined source apportionment method for phosphorus pollution, characterized in that, include: Phosphorus speciation abundance data, driving factor data, and microbial characteristic data consisting of microbial community and phosphorus cycling functional genes were collected from samples from various receptor sampling points and potential pollution sources in the river channel, and phosphorus speciation abundance matrix, driving factor matrix, and microbial characteristic matrix were formed respectively; wherein, the driving factor data includes environmental factor data and relative abundance data of different degradation activity soluble organic matter components. The macroscopic source-sink effect coefficients of each phosphorus form were calculated using the aforementioned driving factor matrix; The biochemical effect coefficients of each phosphorus form were calculated using the aforementioned microbial characteristic matrix; The phosphorus speciation abundance matrix is ​​subjected to a central log-ratio transformation. Then, based on the transformation results, the nonconservative risk and environmental noise ratio of each phosphorus speciation are calculated. Effective tracers are selected according to a preset threshold, and the transformed abundance parameters of the effective tracers are extracted. The calculated macroscopic source-sink effect coefficients of each phosphorus form are coupled with the corresponding biochemical effect coefficients of the same phosphorus form to generate the comprehensive effect coefficients of each phosphorus form. For the selected effective tracers, the initial source spectrum mean and standard deviation of the corresponding pollution sources are synchronously scaled and corrected using the comprehensive effect coefficient, and then re-normalized to obtain the dynamic corrected mean and dynamic corrected standard deviation of each effective tracer, thereby reconstructing the dynamic corrected source spectrum matrix; and the equivalent fractionation coefficient of each pollution source at each sampling point and the corresponding hydrological season is extracted, wherein the equivalent fractionation coefficient is the difference between the dynamic corrected mean and the initial source spectrum mean. The abundance parameters, the dynamically corrected source spectrum matrix, and the equivalent fractionation coefficient are simultaneously input into the Bayesian mixture model. The iterative sampling of the posterior probability distribution is performed through its built-in iterative sampling algorithm. The contribution rate distribution of each pollution source to a specific phosphorus form in the receptor and the corresponding mean are output through mass-weighted calculation. Based on the mean of the posterior probability distribution of the contribution rate of each pollution source and the dynamically corrected source spectrum matrix, the theoretical predicted abundance of each phosphorus species in the receptor is calculated, and the basic residual is obtained by subtracting the theoretical predicted abundance from the actual observed value. A univariate linear regression model is constructed using the theoretically predicted abundance as the independent variable and the actual observed value as the dependent variable. The global analysis accuracy is quantified by the coefficient of determination and root mean square error of the model. The basic residuals are standardized and significance tested to identify significant outliers. If both the global analysis accuracy and the number of significant outliers meet their respective preset thresholds, the source apportionment result is deemed reliable, and the posterior probability distribution of the contribution rate of each pollution source is output. If either the global analysis accuracy or the number of significant outliers does not meet its respective preset threshold, the current source apportionment result is deemed unreliable, and the posterior probability distribution of the contribution rate of each pollution source is not output.

2. The refined source apportionment method for phosphorus pollution as described in claim 1, characterized in that, The phosphorus speciation data includes two main categories: inorganic phosphorus and organic phosphorus. Inorganic phosphorus includes exchangeable phosphorus, aluminum-bound phosphorus, iron-bound phosphorus, and calcium-bound phosphorus. Organic phosphorus includes sodium bicarbonate-extractable organic phosphorus, hydrochloric acid-extractable organic phosphorus, fulvic acid-bound organic phosphorus, humic acid-bound organic phosphorus, and residual organic phosphorus.

3. The refined source apportionment method for phosphorus pollution as described in claim 1, characterized in that, The step of calculating the macroscopic source-sink effect coefficients of each phosphorus form using the driving factor matrix includes: The feature importance score of each driving factor variable in the driving factor matrix to phosphorus form transformation is extracted using the lightweight gradient boosting machine algorithm. The correlation coefficient between each driving factor variable and phosphorus form transformation is calculated by combining the Mantell test. The feature importance score is multiplied by the correlation coefficient to obtain the comprehensive weight of each driving factor variable. The relative abundance of each degradation-active soluble organic matter component is multiplied by the combined weight of the driving factor variable corresponding to that component, and all products are summed to obtain the driving force score. The driving force score is normalized by range normalization to obtain the normalized driving force score. ; Calculate the macroeconomic source-sink effect coefficient based on the normalized driving force score. : ； In the formula: Source-sink direction discrimination coefficient; To correct the step size parameter.

4. The refined source apportionment method for phosphorus pollution as described in claim 1, characterized in that, The step of calculating the biochemical effect coefficients of each phosphorus form using the microbial feature matrix includes: Calculate the global driving weights of microbial community and phosphorus cycling functional genes on each phosphorus form in the microbial feature matrix; The biochemical driving potential of the microbial community and phosphorus cycling functional genes was calculated, and the biochemical effect coefficient of each phosphorus form was obtained by combining the eco-efficiency score.

5. The refined source apportionment method for phosphorus pollution as described in claim 4, characterized in that, Calculate the global driving weights of microbial community and phosphorus cycling functional genes on each phosphorus form in the microbial feature matrix, including: Using the microbial community, phosphorus cycle functional genes, and phosphorus forms as latent variable nodes, the path relationship between the microbial community node and the functional gene node and the phosphorus form node is set. Using a partial least squares path model, the path coefficients from the microbial community nodes to the phosphorus morphology nodes and the path coefficients from the phosphorus cycle functional gene nodes to the phosphorus morphology nodes are extracted, and the sum of the absolute values ​​of the two types of path coefficients is used as the global driving weight.

6. The refined source apportionment method for phosphorus pollution as described in claim 4, characterized in that, The ecological efficiency score was obtained in the following way: A collinear network of phosphorus speciation, microbial community, and functional genes was constructed. The proportion of negatively correlated connections in the network was extracted, the structural stability efficiency was calculated, and then combined with the seasonal metabolic intensity index to obtain the ecological efficiency score.

7. The method for refined source apportionment of phosphorus pollution as described in any one of claims 4 to 6, characterized in that, The calculation of the biochemical driving potential of the microbial community and phosphorus cycling functional genes, combined with the eco-efficiency score, yields the biochemical effect coefficients of each phosphorus form, including: Using the random forest algorithm, with the relative abundance of microbial genera and the relative abundance of phosphorus cycling functional genes in the microbial community as input features and the phosphorus form transformation intensity as the prediction target, the feature contribution weights of each microbial genera and each phosphorus cycling functional gene to phosphorus form transformation are extracted. The correlation coefficient between the relative abundance of the phosphorus cycle functional genes and the phosphorus form conversion intensity was calculated by correlation analysis, and the driving potential of the functional genes was calculated based on the correlation coefficient. The microbial driving potential is calculated based on the correlation coefficient between the relative abundance of the microbial genus and the phosphorus form transformation intensity, as well as the microbial component in the feature contribution weight. The microbial driving potential is multiplied by the functional gene driving potential and scaled to the range of 0.5 to 1.5 using a linear mapping to obtain the standardized comprehensive biochemical driving potential. Multiply the global driving weight, the comprehensive biochemical driving potential, and the ecological efficiency score to obtain the biochemical incremental score. Then, multiply the biochemical incremental score by the biochemical scaling step size parameter and add 1 to obtain the biochemical effect coefficient of each phosphorus form.

8. The refined source apportionment method for phosphorus pollution as described in claim 7, characterized in that, The formula for calculating the functional gene drive potential is as follows: ; In the formula, The functional gene driving potential corresponding to the i-th phosphorus form, denoted as the correlation coefficient between the relative abundance of the j-th phosphorus cycling functional gene and the phosphorus form conversion intensity. Let be the average relative abundance of the j-th phosphorus cycle functional gene in all samples, and m be the total number of phosphorus cycle functional genes.

9. The refined source apportionment method for phosphorus pollution as described in claim 7, characterized in that, The formula for calculating the microbial drive potential is as follows: )； In the formula, The microbial driving potential corresponding to the i-th phosphorus form, is the correlation coefficient between the relative abundance of the k-th microbial genus and the phosphorus form transformation intensity. Let l be the feature contribution weight of the k-th microbial genus to phosphorus form transformation obtained by the random forest algorithm, where l is the total number of microbial genera.

10. A refined source apportionment system for phosphorus pollution, used to implement the refined source apportionment method for phosphorus pollution as described in any one of claims 1 to 9, characterized in that, include: The data acquisition module is configured to collect phosphorus speciation abundance data, driving factor data, and microbial characteristic data composed of microbial community and phosphorus cycle functional genes from samples collected from various receptor sampling points and potential pollution sources in the river channel, forming phosphorus speciation abundance matrix, driving factor matrix, and microbial characteristic matrix, respectively. The macro-source-sink effect quantification module is configured to calculate the macro-source-sink effect coefficients of each phosphorus form using the driving factor matrix. The microscopic biochemical effect quantification module is configured to calculate the biochemical effect coefficients of each phosphorus form using the microbial characteristic matrix; The tracer screening module is configured to perform a central log-ratio transformation on the phosphorus speciation abundance matrix, then calculate the nonconservative risk and environmental noise ratio of each phosphorus speciation based on the transformation result, screen out effective tracers according to a preset threshold, and extract the transformed abundance parameters of the effective tracers. The comprehensive effect coefficient generation module is configured to couple the calculated macroscopic source-sink effect coefficients of each phosphorus form with the corresponding biochemical effect coefficients of the same phosphorus form to generate the comprehensive effect coefficients of each phosphorus form. The source spectrum dynamic correction module is configured to synchronously scale and correct the initial source spectrum mean and standard deviation of the corresponding pollution source for the selected effective tracers using the comprehensive effect coefficient, and then re-normalize them to obtain the dynamic correction mean and dynamic correction standard deviation of each effective tracer, thereby reconstructing the dynamic correction source spectrum matrix. The equivalent fractionation coefficients of each pollution source at each sampling point and in the corresponding hydrological season are extracted. The equivalent fractionation coefficients are the difference between the dynamically corrected mean and the initial source spectrum mean. The Bayesian contribution rate calculation module is configured to simultaneously input the abundance parameters, dynamically corrected source spectrum matrix, and equivalent fractionation coefficients into the Bayesian mixture model, perform iterative sampling of the posterior probability distribution using an iterative sampling algorithm, and output the contribution rate distribution of each pollution source to a specific phosphorus form in the receptor through mass-weighted calculation; and, The result verification module is configured to calculate the theoretical predicted abundance of each phosphorus species in the receptor based on the mean of the posterior probability distribution of the contribution rate of each pollution source and the dynamically corrected source spectrum matrix, and obtain the basic residual by subtracting the theoretical predicted abundance from the actual observed value; a univariate linear regression model is constructed with the theoretical predicted abundance as the independent variable and the actual observed value as the dependent variable, and the global analytical accuracy is quantified by the determination coefficient and root mean square error of the model. The basic residuals are standardized and significance tested to identify significant outliers. If both the global resolution accuracy and the number of significant outliers meet their respective preset thresholds, the source apportionment result is deemed reliable, and the posterior probability distribution of the contribution rate of each pollution source is output. If either the global resolution accuracy or the number of significant outliers does not meet its respective preset threshold, the current source apportionment result is deemed unreliable, and the posterior probability distribution of the contribution rate of each pollution source is not output.

11. A computer-readable storage medium, characterized in that, A computer program is stored on the computer-readable storage medium, which, when executed by a processor, implements the phosphorus pollution refined source apportionment method as described in any one of claims 1 to 9.

Images