Method and system for assessing pollution risk based on environmental microbiome features

Abstract

The application discloses a pollution risk assessment method and system based on environmental microbiome characteristics, and relates to the technical field of environmental pollution risk assessment, and comprises the following steps: obtaining the classification unit abundance and functional gene abundance of a target water body microbial community sample through sequencing, constructing a correlation matrix, and screening a pollution response feature based on the structural difference of adjacent time points; adopting nonlinear mapping to calculate the weighted deviation degree of the feature relative to the environmental quality standard limit value as a risk score; performing time series analysis on the score to obtain an evolution trend value and a period parameter, extracting geometric feature parameters of the evolution track, and dividing a risk grade according to the geometric feature parameters; performing spatial interpolation on the sampling position coordinates to generate a pollution risk spatial distribution map and an evaluation report, so that accurate identification, grading and visualization of environmental pollution risk are realized.

Description

Technical Field

[0001] This invention relates to the field of environmental pollution risk assessment technology, specifically to a pollution risk assessment method and system based on environmental microbiome characteristics. Background Technology

[0002] Traditional pollution risk assessment methods primarily rely on monitoring chemical indicators, such as physicochemical parameters like heavy metal content and organic pollutant concentrations, and assess the degree of pollution by comparing them with environmental quality standards. These methods suffer from limitations such as poor timeliness, low sensitivity, and inability to reflect the true state of the ecosystem. Particularly for complex aquatic ecosystems, a single chemical indicator is insufficient to comprehensively reflect the combined impact of pollutants on the ecosystem and its potential risks.

[0003] Existing microbiome-based environmental assessment methods primarily focus on changes in community structure, lacking in-depth analysis of the correlation between microbial function and environmental pollution, making it difficult to establish a quantitative relationship between microbial characteristics and pollution risk. Furthermore, the analysis of microbial community data is mostly concentrated on static comparisons at a single point in time, lacking spatiotemporal evolution analysis of dynamic changes in microbial communities. When classifying risk levels, existing methods often employ simple linear models or expert judgment, failing to accurately capture the complex nonlinear relationship between environmental pollution and microbial response. The spatial distribution characteristics of pollution risk are also frequently overlooked, resulting in risk assessment results that are insufficient to provide accurate spatial decision support for regional environmental management.

[0004] In screening pollution response characteristics, existing methods mostly rely on differential or correlation analysis, neglecting the intrinsic connections between microbial taxa and functional genes. This simplistic approach fails to fully tap the complex ecological network information contained in microbiome data, reducing the accuracy and predictive power of pollution risk assessment. Summary of the Invention

[0005] The purpose of this invention is to provide a pollution risk assessment method and system based on environmental microbiome characteristics, aiming to solve at least one of the technical problems existing in the prior art.

[0006] The technical solution of this invention is: a pollution risk assessment method based on environmental microbiome characteristics, comprising the following steps: Microbial community samples from the target water body at multiple time points were obtained and sequenced. The sequencing results were then annotated for classification and function to obtain the abundance of microbial taxa and functional genes. The correlation coefficient between the abundance of microbial taxa and the abundance of functional genes was calculated to construct a correlation matrix. Pollution response characteristics were obtained by screening the structural differences in the correlation matrix between adjacent time points. Obtain the environmental quality standard limits corresponding to the pollution response characteristics, and calculate the weighted deviation of the pollution response characteristics from the environmental quality standard limits through nonlinear mapping as the pollution risk score; Time series fitting and periodic analysis are performed on the pollution risk score to obtain the pollution evolution trend value and pollution cycle parameter. Based on the pollution risk score and pollution evolution trend value, the evolution trajectory is constructed and geometric feature parameters are extracted. The pollution risk level is divided according to the geometric feature parameter and pollution cycle parameter. The sampling location coordinates of the microbial community samples are obtained as the interpolation location. Spatial interpolation is performed using the pollution risk score and pollution evolution trend value as the interpolation values. Combined with the pollution risk level, a spatial distribution map of pollution risk is generated and a pollution risk assessment report is output.

[0007] Microbial community samples from the target water body at multiple time points were obtained and sequenced. The sequencing results were then annotated for taxonomic and functional data, yielding the abundance of microbial taxa and functional genes, including: Water samples were collected from the target water body at multiple time points as microbial community samples, and nucleic acids were extracted from the microbial community samples for high-throughput sequencing to obtain sequencing results. Sequence clustering is performed on the sequencing results. Sequences with similarity exceeding a preset similarity threshold are grouped into the same microbial taxonomic unit to complete the classification annotation. The abundance of the initial microbial taxonomic unit is obtained by calculating the ratio of the number of sequences in each microbial taxonomic unit to the total number of sequences in the sequencing results. The frequency of occurrence of each microbial taxonomic unit at multiple time points was statistically analyzed, and the initial microbial taxonomic unit abundance corresponding to the microbial taxonomic units with an occurrence frequency higher than a preset frequency threshold was selected and normalized to obtain the microbial taxonomic unit abundance. The functional gene types encoded by each sequence in the sequencing results are identified. The metabolic pathways to which the functional gene types belong are grouped. The abundance of the metabolic pathway is obtained by summing the proportion of the number of sequences contained in each metabolic pathway to the total number of sequences. The functional gene abundance is obtained by normalization.

[0008] A correlation matrix was constructed by calculating the correlation coefficient between the abundance of microbial taxa and the abundance of functional genes. Pollution response characteristics were then selected based on the structural differences in the correlation matrix between adjacent time points. For each time point, the time series correlation coefficient between the abundance of each microbial taxonomic unit and the abundance of each functional gene is calculated to obtain the correlation coefficient. All correlation coefficients are then filled into a correlation matrix according to the correspondence between microbial taxonomic units and functional genes. Singular value decomposition is performed on the correlation matrix at each time point to obtain the dominant singular values ​​and the corresponding dominant singular vectors. Vector projection is performed on the dominant singular vectors of adjacent time points to calculate the projection residuals. The projection residuals and the changes in the dominant singular values ​​are jointly weighted to obtain the structural difference values ​​of the correlation matrix between adjacent time points, and a sequence of structural difference values ​​is constructed. Baseline fitting is performed on the structural difference value sequence to identify structural anomaly time points that deviate from the baseline. The correlation matrix corresponding to the structural anomaly time point is extracted and the correlation matrix of the baseline time period is calculated to obtain the difference matrix. The difference matrix is ​​divided into blocks and clustered according to the numerical distribution of matrix elements to obtain response modules. The concentration of matrix elements within each response module is calculated, and the microbial taxonomic unit and functional gene corresponding to the response module with the highest concentration are extracted as pollution response features.

[0009] The environmental quality standard limits corresponding to pollution response characteristics are obtained, and the weighted deviation of the pollution response characteristics from the environmental quality standard limits is calculated using a nonlinear mapping to obtain the pollution risk score, which includes: Functional genes are extracted from pollution response characteristics, the metabolic pathways to which the functional genes belong are identified, the integrity of the metabolic pathways is calculated, pollutant types are screened based on the integrity of the metabolic pathways, and the environmental quality standard limits corresponding to the pollutant types are obtained. The abundance of microbial taxa and functional genes corresponding to pollution response characteristics are extracted. The normalized ratios of microbial taxa abundance and functional gene abundance relative to environmental quality standard limits are calculated respectively. The deviation of pollution response characteristics from environmental quality standard limits is obtained by nonlinear mapping transformation of the normalized ratios. Based on the correlation coefficients of pollution response characteristics in the correlation matrix, the contribution rate of the correlation coefficients is calculated to construct a weight vector. The deviation is then weighted and summed using the weight vector to obtain the weighted deviation as the pollution risk score.

[0010] Time-series fitting and periodic analysis are performed on pollution risk scores to obtain pollution evolution trend values ​​and pollution cycle parameters. Based on the pollution risk scores and pollution evolution trend values, an evolution trajectory is constructed and geometric feature parameters are extracted. Pollution risk levels are then classified according to the geometric feature parameters and pollution cycle parameters, including: A time series of pollution risk scores is constructed, the time series is segmented by a sliding window, the average rate of change within each segment is calculated, abrupt changes in the average rate of change between adjacent segments are identified, and the time period is re-divided based on the abrupt changes. The weighted combination of the average rate of change in each time period is calculated as the pollution evolution trend value. Wavelet decomposition is performed on the time series to extract the periodic components, and the dominant oscillation period of the periodic components is calculated as the pollution period parameter. The pollution risk scores are arranged in chronological order as a state vector, and the pollution evolution trend values ​​are arranged in chronological order as a driving vector. The projection trajectory of the state vector in the direction of the driving vector is calculated, and the oscillation amplitude and phase drift of the projection trajectory are extracted as geometric feature parameters through morphological analysis. Risk fluctuation intensity is constructed based on the oscillation amplitude in the geometric characteristic parameters, and periodic stability index is constructed based on the pollution cycle parameters. The coupling index is calculated based on the risk fluctuation intensity and periodic stability index. Multi-level thresholds are set according to the numerical distribution of the coupling index, and the samples are divided into different pollution risk levels according to the multi-level thresholds.

[0011] Calculate the projected trajectory of the state vector in the direction of the driving vector, and perform morphological analysis on the projected trajectory to extract the oscillation amplitude and phase drift of the trajectory, including: Calculate the unit direction vector of the driving vector, calculate the dot product of the state vector and the unit direction vector to obtain the projection scalar sequence, and connect the projection scalar sequence in time order to form the projection trajectory. The upper and lower envelopes are obtained by fitting the projection trajectory with an envelope. The instantaneous distance sequence between the upper and lower envelopes is calculated. The oscillation amplitude is obtained by performing a time-weighted average on the instantaneous distance sequence. The extreme points of the projected trajectory are detected to obtain peak points and valley points. The time interval between adjacent peak points is calculated to obtain the peak interval sequence. The average value of the peak interval sequence is calculated, and the difference between the average value and the pollution cycle parameter is used as the phase drift.

[0012] The sampling location coordinates of the microbial community samples are obtained as the interpolation location. Spatial interpolation is performed using the pollution risk score and pollution evolution trend value as interpolation values. Combined with the pollution risk level, a spatial distribution map of pollution risk is generated, and a pollution risk assessment report is output, including: The sampling location coordinates of the microbial community sample are obtained as the interpolation location, and the interpolation location is divided into spatial grids to obtain grid nodes; Calculate the spatial distance between grid nodes and interpolation locations, construct distance weights based on the spatial distances, use pollution risk scores and pollution evolution trend values ​​as interpolation values, and use the distance weights to perform spatial interpolation on the interpolation values ​​to obtain the interpolated pollution risk scores and interpolated pollution evolution trend values ​​for each grid node. Calculate the difference between the interpolated pollution risk score and the interpolated pollution evolution trend value for each grid node, and identify the spatial location of the risk transition zone based on the difference. Obtain the pollution risk level corresponding to the interpolation location, draw the pollution risk level boundary based on the pollution risk level, and overlay the spatial location of the risk transition zone, the pollution risk level boundary and the interpolated pollution risk score to generate a pollution risk spatial distribution map. The distribution differences of each pollution risk level inside and outside the risk transition zone are statistically analyzed. Based on the distribution differences, the spatial transition patterns of pollution risk levels are identified. The spatial transition patterns, interpolated pollution evolution trend values, and pollution risk spatial distribution maps are integrated to construct and output a pollution risk assessment report.

[0013] This invention provides a pollution risk assessment system based on environmental microbiome characteristics, the system comprising: The sequencing and annotation module is used to obtain microbial community samples of the target water body at multiple time points for sequencing, and to perform classification annotation and functional annotation on the sequencing results to obtain the abundance of microbial taxonomic units and functional genes. The feature screening module is used to calculate the correlation coefficient between the abundance of microbial taxonomic units and the abundance of functional genes to construct a correlation matrix, and to screen pollution response features based on the structural differences in the correlation matrix between adjacent time points. The risk scoring module is used to obtain the environmental quality standard limits corresponding to the pollution response characteristics. It calculates the weighted deviation of the pollution response characteristics from the environmental quality standard limits through nonlinear mapping as the pollution risk score. The classification module is used to perform time-series fitting and periodic analysis on pollution risk scores to obtain pollution evolution trend values ​​and pollution cycle parameters. Based on the pollution risk scores and pollution evolution trend values, the module constructs evolution trajectories and extracts geometric feature parameters. Based on the geometric feature parameters and pollution cycle parameters, the module classifies pollution risk levels. The spatial analysis module is used to obtain the sampling location coordinates of microbial community samples as interpolation locations, use pollution risk scores and pollution evolution trend values ​​as interpolation values ​​for spatial interpolation, combine pollution risk levels to generate a pollution risk spatial distribution map and output a pollution risk assessment report.

[0014] One technical solution provided in this embodiment of the invention is an electronic device, including: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the steps in any of the aforementioned methods.

[0015] One technical solution provided in this embodiment of the invention is a computer-readable storage medium storing computer program instructions, which, when executed by a processor, implement the steps in any of the aforementioned methods.

[0016] This invention captures early warning signals of environmental pollution by analyzing the temporal dynamic changes of microbiomes, thereby improving the sensitivity and timeliness of pollution risk identification. By combining the correlation analysis of microbial taxonomic units and functional genes, it screens out highly responsive pollution characteristic indicators, enhancing the specificity of risk assessment. Using a nonlinear mapping method, it establishes a quantitative relationship between microbial characteristics and environmental quality standards, achieving quantitative scoring of pollution risks. Through temporal fitting and periodic analysis, it reveals the dynamic patterns and evolutionary trends of pollution changes, improving the accuracy of risk prediction. The risk level classification method based on geometric characteristic parameters and pollution cycle parameters makes the assessment results more objective and scientific, realizing spatial visualization of pollution risks and providing precise spatial decision support for regional environmental management. Attached Figure Description

[0017] Figure 1 A flowchart illustrating a pollution risk assessment method based on environmental microbiome characteristics provided in an embodiment of the present invention; Figure 2 This is a schematic diagram of the pollution risk assessment system based on environmental microbiome characteristics according to an embodiment of the present invention. Detailed Implementation

[0018] like Figure 1 As shown, Figure 1 A flowchart of a pollution risk assessment method based on environmental microbiome characteristics provided in an embodiment of the present invention is shown. The method includes the following steps: Step 101: Obtain microbial community samples from the target water body at multiple time points for sequencing, and perform classification annotation and functional annotation on the sequencing results to obtain the abundance of microbial taxonomic units and functional genes.

[0019] In some embodiments of the present invention, step 101 may specifically include the following sub-steps: Sub-step 1011: Collect water samples from the target water body at multiple time points as microbial community samples, extract nucleic acids from the microbial community samples and perform high-throughput sequencing to obtain sequencing results; Sub-step 1012: Perform sequence clustering on the sequencing results, group sequences with sequence similarity exceeding a preset similarity threshold into the same microbial taxonomic unit to complete the classification annotation, and calculate the ratio of the number of sequences in each microbial taxonomic unit to the total number of sequences in the sequencing results to obtain the initial microbial taxonomic unit abundance. Sub-step 1013: Count the frequency of occurrence of each microbial taxonomic unit at multiple time points, screen the initial microbial taxonomic unit abundance corresponding to the microbial taxonomic units whose occurrence frequency is higher than the preset frequency threshold, and normalize the microbial taxonomic unit abundance. Sub-step 1014: Identify the functional gene types encoded by each sequence in the sequencing results, group them according to the metabolic pathways to which the functional gene types belong, sum the ratio of the number of sequences contained in each metabolic pathway to the total number of sequences to obtain the metabolic pathway abundance, and normalize the results to obtain the functional gene abundance.

[0020] First, microbial community samples from the target water body at different time points are obtained and sequenced. In practice, water samples are collected from the target water body as microbial community samples. Typically, sterile sampling bottles are used to collect samples at a depth of 20 cm below the water surface, with at least 500 mL of sample collected at each sampling point. The collected water samples are filtered through a filter membrane, and the microbial cells on the filter membrane are collected. A nucleic acid extraction kit is used to extract total nucleic acids, including DNA and RNA, from the microbial community samples. The extracted nucleic acids need to be purified to remove impurities that may affect subsequent sequencing. After the purified nucleic acid samples pass quality testing, sequencing libraries are constructed, followed by high-throughput sequencing. Currently commonly used sequencing technologies include amplicon sequencing and metagenomic sequencing. The former mainly targets specific marker genes such as 16S rRNA genes, while the latter sequences the entire genome in the sample.

[0021] The obtained sequencing data undergoes quality control to remove low-quality and chimeric sequences. Sequences that pass quality control are then subjected to cluster analysis, grouping sequences with similarity exceeding a preset threshold into the same microbial taxonomic unit. Typically, 97% similarity is used as the threshold for species-level classification. Sequence clustering can be achieved using distance-based or seed-based clustering methods. After clustering, the central sequences of each cluster are compared with sequences in a reference database to determine the taxonomic information represented by each cluster, providing hierarchical annotation from phylum, class, order, family, genus to species.

[0022] Based on clustering and annotation results, the proportion of sequences contained in each microbial taxonomic unit to the total number of sequences was calculated to obtain the initial abundance of the microbial taxonomic units. To eliminate the interference of incidentally occurring microbial taxonomic units on the analysis results, the frequency of occurrence of each microbial taxonomic unit in samples at different time points was calculated. An appropriate frequency threshold was set; for example, microbial taxonomic units appearing in more than 50% of the samples were considered core microbial taxonomic units. The initial abundance data corresponding to the core microbial taxonomic units were screened out, and these data were normalized. The normalization process can use the relative abundance transformation method, with the relative abundance value x = a / m·d, where a represents the number of sequences in each microbial taxonomic unit, m represents the total number of sequences in the sample, and d is a fixed coefficient.

[0023] In addition to classification annotation, functional annotation of sequencing results is also necessary to identify the functional gene types encoded by each sequence. Functional annotation can be accomplished by comparing sequences with functional gene databases, commonly used functional gene databases include protein function databases and metabolic pathway databases. Sequence homology searches determine the function of the protein encoded by each sequence and its location within a metabolic pathway. Functional genes are then grouped according to their corresponding metabolic pathways, such as carbon cycle, nitrogen cycle, sulfur cycle, phosphorus cycle, etc. For each metabolic pathway, the abundance of the pathway is calculated by summing the proportions of all functional gene sequences contained within that pathway relative to the total number of sequences.

[0024] The abundance of metabolic pathways is normalized to obtain the final functional gene abundance. Normalization can be performed using logarithmic transformation or standardization, with the aim of reducing data range differences and making the abundance data of different metabolic pathways comparable. The normalized functional gene abundance reflects the potential of the microbial community to perform various biogeochemical processes and can be used for subsequent pollution risk assessment analysis.

[0025] This invention establishes a pollution risk assessment technology system based on environmental microbiome characteristics through in-depth analysis of the target aquatic microbial community. Compared with traditional physicochemical indicator assessment methods, this method can detect environmental pollution risks earlier and more sensitively, reflecting changes in ecosystem health. By accurately measuring the abundance of microbial taxa and functional genes, it is possible not only to monitor the presence of pollutants but also to assess their actual impact on ecosystem function.

[0026] Step 102: Calculate the correlation coefficient between the abundance of microbial taxonomic units and the abundance of functional genes to construct a correlation matrix. Based on the structural differences in the correlation matrix between adjacent time points, pollution response characteristics are screened.

[0027] In some embodiments of the present invention, step 102 may specifically include the following sub-steps: Sub-step 1021: For each time point, calculate the time series correlation coefficient between the abundance of each microbial taxonomic unit and the abundance of each functional gene, and fill all the correlation coefficients into an association matrix according to the correspondence between microbial taxonomic units and functional genes. Sub-step 1022: Perform singular value decomposition on the correlation matrix at each time point to obtain the dominant singular value and the corresponding dominant singular vector. Perform vector projection on the dominant singular vectors of adjacent time points to calculate the projection residual. Perform joint weighting on the projection residual and the change of dominant singular value to obtain the structural difference value of the correlation matrix between adjacent time points, and construct a sequence of structural difference values. Sub-step 1023: Baseline fitting is performed on the structural difference value sequence to identify structural anomaly time points that deviate from the baseline. The correlation matrix corresponding to the structural anomaly time point is extracted and the correlation matrix of the baseline time period is calculated to obtain the difference matrix. Sub-step 1024: The difference matrix is ​​divided into blocks and clustered according to the numerical distribution of matrix elements to obtain response modules. The concentration of matrix elements within each response module is calculated, and the microbial taxonomic unit and functional gene corresponding to the response module with the highest concentration are extracted as pollution response features.

[0028] For samples collected at each time point, the time-series correlation between the abundance of each microbial taxonomic unit and the abundance of each functional gene was calculated. The time-series correlation could be calculated using Pearson correlation coefficient or Spearman rank correlation coefficient. Pearson correlation coefficient is suitable for linear relationships, while Spearman rank correlation coefficient is suitable for non-linear relationships. The correlation coefficient between each pair of microbial taxonomic units and functional genes was calculated, with values ​​ranging from -1 to 1. Values ​​close to 1 indicate a strong positive correlation, values ​​close to -1 indicate a strong negative correlation, and values ​​close to 0 indicate no significant correlation.

[0029] The calculated association coefficients are filled into an association matrix according to the correspondence between microbial taxonomic units and functional genes. The rows of this matrix represent microbial taxonomic units, the columns represent functional genes, and each element represents the association coefficient between the corresponding taxonomic unit and functional gene. For easier subsequent analysis, the association matrix can be visualized, such as using a heatmap to show the differences in association strength.

[0030] Singular value decomposition (SVD) is performed on the correlation matrix constructed at each time point, decomposing the matrix into a product of three matrices: the left singular vector matrix, the singular value diagonal matrix, and the transpose of the right singular vector matrix. The singular values ​​are arranged in descending order, with the largest singular value called the dominant singular value, and its corresponding left and right singular vectors called the dominant singular vectors. The dominant singular value and its corresponding dominant singular vectors contain the main structural features of the correlation matrix.

[0031] The vector projection between the dominant singular vectors at adjacent time points is calculated. The projection coefficients are obtained by projecting the current dominant singular vector onto the previous dominant singular vector. The projection residual is obtained by subtracting the projected portion from the current dominant singular vector. The magnitude of the projection residual reflects the degree of change in the direction of the dominant singular vectors at adjacent time points. Simultaneously, the change in the dominant singular value at adjacent time points is calculated, i.e., the difference between the current dominant singular value and the previous dominant singular value. The projection residual and the change in the dominant singular value are jointly weighted to obtain the structural difference value of the incidence matrix between adjacent time points. The joint weighting can use a linear combination of the projection residual and the change in singular value, and the weights can be adjusted according to the actual application scenario.

[0032] Structural discrepancies between adjacent time points are arranged chronologically to form a structural discrepancy sequence. Baseline fitting is then performed on this sequence to identify structurally anomalous time points that significantly deviate from the baseline. Baseline fitting can employ techniques such as moving averages or multinomial regression. An appropriate threshold is determined; when a structural discrepancy exceeds a certain multiple of the baseline standard deviation, that time point is identified as a structurally anomalous.

[0033] The correlation matrix corresponding to the structural anomaly time point is compared with the correlation matrix of the baseline time period. The difference matrix is ​​obtained by calculating the difference between the matrices. The baseline time period can be a continuous period with stable and low structural difference values. The average value of each correlation matrix within this period is calculated as the baseline correlation matrix. Each element of the difference matrix represents the magnitude of change in the correlation strength between the corresponding microbial taxonomic unit and functional gene at the structural anomaly time point compared to the baseline time period.

[0034] The difference matrix is ​​divided into blocks and clustered according to the numerical distribution of its elements to identify response modules. Spectral clustering algorithms can be used for this block clustering, grouping matrix elements with similar values ​​and variation patterns into the same response module. The concentration of the matrix elements within each response module is calculated. Concentration can be defined as the ratio of the standard deviation of elements within the module to the standard deviation of all elements in the matrix, or the ratio of the mean of elements within the module to the mean of all elements in the matrix.

[0035] The response module with the highest concentration was selected, and the corresponding microbial taxa and functional genes were extracted as pollution response features. These features represent microbial groups and functional pathways that show significant correlational changes when environmental pollution occurs, and can serve as indicators for assessing environmental pollution risk. To verify the reliability of the extracted features, the consistency of pollution response features obtained at different structural anomaly time points can be compared; features with high consistency tend to have better indicative value.

[0036] This invention achieves accurate assessment of environmental pollution risks by analyzing the dynamic changes in the network of relationships between microbial community structure and function. Compared to methods based solely on a single microbial community or functional gene, network analysis can capture the overall synergistic patterns of microbial community responses to environmental changes, providing more comprehensive information on ecosystem health. It can extract the most sensitive indicators of pollution from complex microbiome data, reducing the complexity and cost of environmental monitoring and improving the sensitivity and specificity of pollution risk assessment. By monitoring changes in these pollution response characteristics in real time, early warning of environmental pollution can be achieved, providing a scientific basis for environmental management and ecological protection.

[0037] Step 103: Obtain the environmental quality standard limit corresponding to the pollution response characteristics, and calculate the weighted deviation of the pollution response characteristics from the environmental quality standard limit as the pollution risk score through nonlinear mapping.

[0038] In some embodiments of the present invention, step 103 may specifically include the following sub-steps: Sub-step 1031: Extract functional genes from pollution response characteristics, identify the metabolic pathways to which the functional genes belong, calculate the integrity of the metabolic pathways, screen pollutant types based on the integrity of the metabolic pathways, and obtain the environmental quality standard limits corresponding to the pollutant types. Sub-step 1032: Extract the abundance of microbial taxonomic units and functional gene abundance corresponding to the pollution response characteristics, calculate the normalized ratio of microbial taxonomic unit abundance relative to the environmental quality standard limit and the normalized ratio of functional gene abundance relative to the environmental quality standard limit, and perform nonlinear mapping transformation on the normalized ratio to obtain the deviation of the pollution response characteristics from the environmental quality standard limit. Sub-step 1033: Based on the correlation coefficients of the pollution response characteristics in the correlation matrix, calculate the contribution rate of the correlation coefficients to construct a weight vector, and use the weight vector to perform a weighted summation of the deviations to obtain the weighted deviation as the pollution risk score.

[0039] The pollution response characteristics obtained through the aforementioned steps include microbial taxonomic units and functional genes. The functional genes need to be isolated for further analysis. The extracted functional genes are annotated by querying functional gene databases to identify the metabolic pathway to which each gene belongs. Common metabolic pathways include biochemical processes related to the metabolism of environmental pollutants, such as the carbon cycle, nitrogen cycle, sulfur cycle, aromatic hydrocarbon degradation, and heavy metal transformation.

[0040] Metabolic pathway integrity is an indicator of the extent to which a set of functional genes covers a specific metabolic process. It is calculated by dividing the number of key enzyme-coding genes detected in a target metabolic pathway by the total number of key enzyme-coding genes required for that pathway. By calculating the integrity of each metabolic pathway, the types of pollutants that may be present in the environment can be inferred. High metabolic pathway integrity means that the microbial community has the complete ability to metabolize that type of pollutant, indicating the potential presence of that pollutant in the environment. A threshold is set based on metabolic pathway integrity to screen out the pollutant types corresponding to metabolic pathways with integrity exceeding the threshold.

[0041] For the pollutant types identified through screening, consult the environmental quality standard limits. Environmental quality standard limits refer to the maximum allowable concentration of a specific pollutant in the environment, typically derived from relevant environmental protection standards. Different types of water bodies, such as surface water, groundwater, and drinking water sources, have different standard limit requirements. When obtaining standard limits, the intended use and functional zoning of the target water body must be considered, and applicable international standard limits should be selected.

[0042] The abundance data of corresponding microbial taxa and functional genes were extracted from the pollution response characteristics. The normalized ratios of microbial taxa abundance relative to environmental quality standard limits and the normalized ratios of functional gene abundance relative to environmental quality standard limits were calculated, respectively. When calculating the normalized ratios, the ratios of abundance values ​​to standard limits were standardized to a uniform range, ensuring comparability between different indicators.

[0043] A nonlinear mapping transformation is applied to the normalized ratio to obtain the deviation of the pollution response characteristics from the environmental quality standard limits. The nonlinear mapping can employ exponential functions or hyperbolic tangent functions to convert the normalized ratio into a deviation. The purpose of the nonlinear mapping is to emphasize the portion exceeding the standard limits, making the risk assessment more sensitive to standard exceedances. The pollutant standard limits corresponding to the target water body type are obtained, and the pollutants are classified into highly toxic and low-toxic pollutants. The choice of mapping function should consider the hazardous characteristics of the pollutants; a steep nonlinear function can be used for highly toxic pollutants, while a gentler nonlinear function can be used for low-toxic pollutants.

[0044] A weight vector is calculated based on the correlation coefficients of pollution response features in the correlation matrix. The correlation coefficients reflect the strength of the interaction between microbial taxa and functional genes; a higher correlation coefficient indicates a stronger interaction and a greater contribution to pollution risk assessment. To calculate the contribution rate of the correlation coefficients, the absolute value of each coefficient is divided by the sum of the absolute values ​​of all correlation coefficients to obtain a normalized contribution rate. A weight vector is then constructed based on the contribution rate, with each element in the weight vector corresponding to a weight value for a pollution response feature.

[0045] The deviations are weighted and summed using a weight vector to obtain a weighted deviation score, which serves as the pollution risk score. The dot product of the weight vector and the deviation vector is the weighted deviation score. This weighting operation ensures that features contributing more to the pollution risk assessment have higher weights, thus improving the accuracy of the risk assessment. The resulting pollution risk score is a comprehensive indicator that reflects the degree of impact of pollutants on the microbial community and its deviation from environmental quality standards.

[0046] Pollution risk scores can be categorized into different risk levels, such as low risk, medium risk, and high risk. The classification of risk levels should consider historical data and expert experience, and set reasonable thresholds. Risk levels can guide environmental management departments to take corresponding prevention and control measures; low-risk areas can continue monitoring, medium-risk areas require increased monitoring frequency, and high-risk areas require immediate pollution control measures.

[0047] This invention achieves precise assessment of water pollution risks by analyzing the response characteristics of the environmental microbiome to pollution and combining this with environmental quality standard limits. Compared to traditional physicochemical indicator monitoring methods, it can detect potential risks before pollutant concentrations reach standard limits, providing early warning for environmental management. This invention integrates microbial ecology, environmental science, and data analysis techniques, establishing a technical pathway from microbiome characteristics to quantitative assessment of pollution risks. It can be applied not only to conventional pollutant monitoring but also to risk assessment of emerging and mixed pollutants, improving the comprehensiveness and foresight of environmental monitoring.

[0048] Step 104: Perform time series fitting and periodic analysis on the pollution risk score to obtain the pollution evolution trend value and pollution cycle parameter. Construct the evolution trajectory based on the pollution risk score and pollution evolution trend value and extract geometric feature parameters. Divide the pollution risk level according to the geometric feature parameter and pollution cycle parameter.

[0049] In some embodiments of the present invention, step 104 may specifically include the following sub-steps: Sub-step 1041: Construct a time series of pollution risk scores, divide the time series into sliding window segments, calculate the average rate of change within each segment, identify abrupt changes in the average rate of change between adjacent segments, and re-divide the time period based on the abrupt changes. Sub-step 1042: Calculate the weighted combination of the average rate of change in each time period as the pollution evolution trend value, perform wavelet decomposition on the time series to extract periodic components, and calculate the dominant oscillation period of the periodic components as the pollution period parameter. Sub-step 1043: Arrange the pollution risk scores in chronological order as a state vector, arrange the pollution evolution trend values ​​in chronological order as a driving vector, calculate the projection trajectory of the state vector in the direction of the driving vector, and perform morphological analysis on the projection trajectory to extract the oscillation amplitude and phase drift of the trajectory as geometric feature parameters. Sub-step 1044: Construct the risk fluctuation intensity based on the oscillation amplitude in the geometric feature parameters, and construct the cycle stability index based on the pollution cycle parameters; Sub-step 1045: Calculate the coupling index based on the risk fluctuation intensity and cycle stability index, set multi-level thresholds according to the numerical distribution of the coupling index, and divide the sample into different pollution risk levels according to the multi-level thresholds.

[0050] Pollution risk scores are arranged chronologically to construct a time series, recording the risk assessment results corresponding to different sampling time points. Sampling time points can be collected at equal intervals, or different sampling frequencies can be set according to actual needs. The constructed time series is segmented using a sliding window. The window size can be adjusted according to the data length and analysis requirements, typically choosing a time scale that can encompass the complete pollution evolution process. The step size of the sliding window can also be flexibly set; a smaller step size can capture more subtle changes but increases computational load; a larger step size is computationally efficient but may ignore short-term fluctuations.

[0051] Within each sliding window, the average rate of change of the pollution risk score is calculated, which is the ratio of the change in risk score within the window to the time interval. If the pollution risk score shows an overall upward trend within the window, the rate of change is positive; if it shows a downward trend, the rate of change is negative. After calculating the average rate of change for all windows, points where the rate of change changes significantly is identified by comparing the differences between adjacent windows; these points are defined as abrupt change points. Abrupt change points can be determined by using a method where the difference in the rate of change exceeds a preset threshold, which can be determined based on the fluctuation range of historical data.

[0052] Based on the identified abrupt change points, the entire time series was redivided into several time periods. Within each time period, the rate of change of pollution risk score was relatively stable, while the rate of change differed significantly between adjacent time periods. This division method can reflect the stage-specific characteristics of the pollution evolution process, such as different stages like pollutant accumulation, stabilization, and decline.

[0053] The pollution evolution trend value is calculated by weighting and combining the average rates of change within the redefined time periods. The weighting can be determined based on the length of the time period or the degree of fluctuation in the pollution risk score within that period; longer time periods or periods with greater fluctuations can be assigned higher weights. The pollution evolution trend value reflects the overall direction and speed of change in pollution risk throughout the entire observation period. Positive values ​​indicate an overall increase in risk, while negative values ​​indicate an overall decrease in risk. The absolute value indicates the rate of change.

[0054] Simultaneously, wavelet decomposition is performed on the pollution risk score time series, breaking it down into wavelets of different frequencies and extracting the periodic components. Morlet wavelets can be used for wavelet decomposition. By analyzing the wavelet energy spectrum, the frequency component with the highest energy is determined as the dominant oscillation period. This period reflects the timescale characteristics of pollution risk changes and can be used as a pollution period parameter.

[0055] Pollution risk scores are arranged chronologically to form a state vector, which describes the pollution risk state of the environmental system at different times. Similarly, pollution evolution trend values ​​are arranged chronologically to form a driving vector, which represents the force driving changes in the environmental system's state. The projection of the state vector onto the driving vector direction yields the projection trajectory. This projection trajectory reflects the dynamic change of the pollution risk state driven by the evolution trend and can be used to analyze the evolutionary patterns of pollution risk.

[0056] Morphological analysis is performed on the projected trajectory to extract its geometric characteristic parameters. Key geometric characteristic parameters include oscillation amplitude and phase drift. Oscillation amplitude refers to the maximum deviation of the projected trajectory from its centerline, reflecting the severity of pollution risk fluctuations. Phase drift refers to the change in phase angle between adjacent peaks or troughs in the projected trajectory, reflecting the stability of the pollution cycle. A larger oscillation amplitude indicates more severe pollution risk fluctuations; a larger phase drift indicates a more unstable pollution cycle.

[0057] A risk fluctuation intensity index is constructed based on the oscillation amplitude in the geometric characteristic parameters. Risk fluctuation intensity can be defined as the ratio of oscillation amplitude to the average oscillation value. This index reflects the relative intensity of pollution risk fluctuations. High risk fluctuation intensity indicates that the environmental system is sensitive to pollution and prone to drastic fluctuations; low risk fluctuation intensity indicates that the environmental system has a certain buffering capacity, and pollution risk changes are relatively gradual.

[0058] A cycle stability index is constructed based on pollution cycle parameters. Cycle stability can be characterized by the coefficient of variation or the standard deviation of the phase drift of the pollution cycle. High cycle stability indicates that the pollution risk changes regularly and is highly predictable; low cycle stability indicates that the pollution risk changes without obvious regularity and has high uncertainty.

[0059] A coupled index is calculated based on the intensity of risk fluctuations and the periodic stability index. This coupled index can be the geometric mean or weighted harmonic mean of the two indices. It comprehensively considers both the fluctuation and periodic characteristics of pollution risk, providing a more complete reflection of its dynamic behavior. Multiple threshold levels are set based on the numerical distribution of the coupled index, using quantile methods or cluster analysis. The samples are then divided into different pollution risk levels, such as low risk, medium risk, high risk, and extremely high risk, according to these threshold levels.

[0060] Different risk levels correspond to different management measures and response strategies. Low risk levels indicate a relatively small risk of environmental pollution, and routine monitoring can be maintained; medium risk levels require increased monitoring frequency and attention to risk trends; high risk levels require the activation of early warning mechanisms and the implementation of proactive prevention and control measures; and extremely high risk levels require immediate emergency response measures to prevent the pollution incident from escalating.

[0061] This invention achieves a comprehensive characterization and precise classification of the dynamic characteristics of environmental pollution risks by performing time-series fitting and periodic analysis on pollution risk scores, combined with the geometric feature parameters of the evolution trajectory. Compared with traditional static risk assessment methods, it can capture the temporal evolution patterns and periodic fluctuation characteristics of pollution risks, improving the dynamic adaptability and predictive ability of risk assessment. It constructs a complete technical path from time-series data to risk levels, enabling the assessment of current pollution risk levels and the prediction of future risk trends.

[0062] Sub-step 1043, calculating the projected trajectory of the state vector in the direction of the driving vector, and performing morphological analysis on the projected trajectory to extract the oscillation amplitude and phase drift of the trajectory, further includes: Calculate the unit direction vector of the driving vector, calculate the dot product of the state vector and the unit direction vector to obtain the projection scalar sequence, and connect the projection scalar sequence in time order to form the projection trajectory. The upper and lower envelopes are obtained by fitting the projection trajectory with an envelope. The instantaneous distance sequence between the upper and lower envelopes is calculated. The oscillation amplitude is obtained by performing a time-weighted average on the instantaneous distance sequence. The extreme points of the projected trajectory are detected to obtain peak points and valley points. The time interval between adjacent peak points is calculated to obtain the peak interval sequence. The average value of the peak interval sequence is calculated, and the difference between the average value and the pollution cycle parameter is used as the phase drift.

[0063] The calculation of the projection trajectory begins with the normalization of the driving vector, converting it into a unit direction vector. This driving vector consists of pollution evolution trend values ​​arranged chronologically, reflecting the dynamic characteristics of pollution risk changes over time. Normalization involves dividing the driving vector by its magnitude, ensuring the resulting unit direction vector has a length of 1 while maintaining its original direction. This process ensures standardization in subsequent projection calculations, facilitating comparisons between different samples. The unit direction vector represents the dominant direction of pollution evolution and serves as an important reference coordinate system for understanding the dynamic evolution of pollution.

[0064] Let the state vector be s = [s1, s2, ..., sn]. n The unit direction vector is u = [u1, u2, ..., u]. n Then their dot product formula is: , where s i It is the i-th element of the state vector, u i It is the i-th element of the unit direction vector.

[0065] The state vector consists of pollution risk scores arranged in chronological order, representing the time-series state of environmental pollution risk. The physical meaning of the dot product operation is to calculate the projection length of the state vector onto the driving vector direction, reflecting the degree to which the pollution risk state is affected by the driving factors. The dot product result forms a sequence of projected scalars, where each scalar value represents the projection intensity of the state onto the driving direction at the corresponding time.

[0066] Projected scalar sequences are connected sequentially over time to form a projected trajectory, maintaining the correspondence between time points during the connection process to ensure the trajectory reflects the true temporal evolution. The projected trajectory can be viewed as the movement trajectory of a pollution risk state under the influence of driving forces, containing key information about the dynamic changes in pollution risk. The shape characteristics of the projected trajectory, such as oscillation frequency, amplitude, and waveform symmetry, can reflect the response pattern and sensitivity of the environmental system to pollution.

[0067] Envelope fitting is performed on the projected trajectory to obtain the upper and lower envelopes. Envelope fitting employs spline interpolation or polynomial fitting methods. The upper envelope is obtained by connecting the local maxima of the trajectory, and the lower envelope is obtained by connecting the local minima. The influence of data noise should be considered during the fitting process. If necessary, the projected trajectory can be smoothed first, such as using moving averages or wavelet denoising techniques, to reduce noise interference with the envelope fitting. The upper and lower envelopes together constitute the oscillation boundary of the projected trajectory, reflecting the range of trajectory fluctuations.

[0068] The instantaneous distance sequence between the upper and lower envelopes is calculated. Instantaneous distance refers to the absolute value of the difference between the upper and lower envelope values ​​at the same time point. The instantaneous distance sequence reflects the change in the oscillation intensity of the projected trajectory at different time points; regions with large instantaneous distances correspond to periods of large trajectory oscillation amplitudes, while regions with small instantaneous distances correspond to periods of small trajectory oscillation amplitudes. The fluctuation characteristics of the instantaneous distance sequence can reveal the time-varying nature of pollution risk fluctuations.

[0069] The oscillation amplitude is obtained by performing a time-weighted average on the instantaneous distance series. This time-weighted averaging considers the differences in importance at different time points, typically assigning higher weights to recent data and lower weights to older data. The weighting function can be an exponential decay function or a linear decay function, with the weights normalized to 1. The oscillation amplitude obtained after time-weighted averaging is a comprehensive indicator, representing the weighted average of the oscillation intensity of the projected trajectory over the entire observation period, reflecting the overall intensity of pollution risk fluctuations.

[0070] Extreme point detection is performed on the projected trajectory. This can be achieved using the derivative test or window comparison method to identify local maxima as peak points and local minima as trough points. To improve detection accuracy, amplitude and duration thresholds can be set to filter out false extreme points caused by noise. Peak and trough points are key characteristic points of the projected trajectory, marking the turning points in the pollution risk state.

[0071] The peak interval sequence is obtained by calculating the time interval between adjacent peak points. The time interval calculation is based on the timestamps corresponding to the peak points, reflecting the periodicity of the projected trajectory oscillation. The statistical properties of the peak interval sequence, such as mean, variance, and coefficient of variation, can be used to analyze the stability and regularity of pollution risk cycles. The changing trend of the peak interval can also reflect the evolution of the pollution cycle; for example, a shortening cycle may indicate an increased pollution risk, while a lengthening cycle may indicate a reduction in pollution.

[0072] The average value of the peak interval sequence is calculated as the average oscillation period of the projected trajectory. The average oscillation period represents the average time required for the projected trajectory to complete one full oscillation, and is a quantitative characterization of the trajectory's periodicity. The average value can be calculated using either the arithmetic mean or the median; the former is sensitive to extreme values, while the latter is robust. For cases with a clear periodic trend, a piecewise averaging method can be used to calculate the average period for different stages.

[0073] The difference between the average oscillation period and the pollution period parameter is used as the phase drift. The pollution period parameter is the dominant oscillation period of the periodic component extracted by wavelet decomposition, representing the inherent periodicity of the pollution risk time series. The phase drift reflects the degree of deviation between the actual oscillation period and the theoretical period of the projected trajectory; a positive value indicates a longer actual period, and a negative value indicates a shorter actual period. The larger the absolute value of the phase drift, the more unstable the periodicity of the pollution risk, and the greater the difficulty in prediction; a phase drift close to zero indicates that the pollution risk has a stable periodicity, which is beneficial for prediction and management.

[0074] The oscillation amplitude and phase drift extracted using the above methods are important geometric characteristic parameters of the dynamic properties of pollution risk, providing a quantitative basis for subsequent risk level classification and management strategy formulation. These parameters can characterize the fluctuation intensity and periodic stability of pollution risk over time, contributing to a deeper understanding of the dynamic evolution of pollution risk.

[0075] This invention achieves dimensionality reduction and feature extraction of complex dynamic processes of pollution risk by projecting multidimensional pollution risk states onto the driving vector direction. This method effectively captures the oscillating characteristics and periodic variation patterns of pollution risk, overcoming the limitations of traditional methods in characterizing dynamic features. The geometric feature parameters extracted by this method intuitively reflect the dynamic change patterns of pollution risk, facilitating understanding and application by environmental management personnel. Through quantitative analysis of oscillation amplitude and phase drift, it achieves accurate assessment of the intensity and periodic stability of pollution risk fluctuations, providing a scientific basis for risk classification.

[0076] Step 105: Obtain the sampling location coordinates of the microbial community sample as the interpolation location, use the pollution risk score and pollution evolution trend value as the interpolation value for spatial interpolation, combine the pollution risk level to generate a pollution risk spatial distribution map and output a pollution risk assessment report.

[0077] In some embodiments of the present invention, step 105 may specifically include the following sub-steps: Sub-step 1051: Obtain the sampling location coordinates of the microbial community sample as the interpolation location, and divide the interpolation location into a spatial grid to obtain grid nodes; Sub-step 1052: Calculate the spatial distance between the grid node and the interpolation location, construct the distance weight based on the spatial distance, use the pollution risk score and pollution evolution trend value as the interpolation value, and use the distance weight to perform spatial interpolation on the interpolation value to obtain the interpolated pollution risk score and interpolated pollution evolution trend value of each grid node. Sub-step 1053: Calculate the difference between the interpolated pollution risk score and the interpolated pollution evolution trend value of each grid node, and identify the spatial location of the risk transition zone based on the difference. Sub-step 1054: Obtain the pollution risk level corresponding to the interpolation position, draw the pollution risk level boundary according to the pollution risk level, and overlay the spatial location of the risk transition zone, the pollution risk level boundary and the interpolated pollution risk score to generate a pollution risk spatial distribution map. Sub-step 1055: Statistically analyze the distribution differences of each pollution risk level inside and outside the risk transition zone, identify the spatial transition patterns of pollution risk levels based on the distribution differences, and integrate the spatial transition patterns, interpolated pollution evolution trend values, and pollution risk spatial distribution maps to construct and output a pollution risk assessment report.

[0078] The sampling location coordinates include longitude, latitude, and elevation information, recording the geographical location of the environmental sample collection. Spatial grid division can employ either equidistant or adaptive strategies. Equidistant division uses the same grid spacing across the entire area; adaptive division adjusts the grid density based on regional characteristics and sample density. Grid nodes serve as interpolation target points, and their location coordinates are determined by the grid division rules.

[0079] The spatial distance between grid nodes and interpolation locations can be calculated using Euclidean distance or a weighted distance that considers terrain factors. Distance weights are constructed based on the inverse relationship between spatial distances, typically employing inverse distance weighting or a Gaussian decay function. Inverse distance weighting uses the reciprocal of the distance as the basis for the weights; the Gaussian decay function controls the weight decay rate by adjusting the bandwidth parameter. Weight normalization ensures that the sum of all weights is 1, maintaining the unbiasedness of the interpolation result.

[0080] Spatial interpolation is performed using pollution risk scores and pollution evolution trend values ​​as interpolation values. Interpolation methods can include inverse distance weighted interpolation, kriging interpolation, or natural neighbor interpolation. Inverse distance weighted interpolation is computationally simple; kriging interpolation considers the autocorrelation of spatial variables; and natural neighbor interpolation is suitable for irregularly distributed samples. During the interpolation process, the influence of environmental factors such as topography and water system distribution on pollution diffusion can be considered, and environmental relevance weights can be introduced to adjust the interpolation results.

[0081] The degree of difference between the interpolated pollution risk score and the interpolated pollution evolution trend value for each grid node is calculated. The degree of difference can be calculated using normalized difference or relative rate of change, reflecting the consistency between the current state of pollution and its future evolution direction. The spatial distribution of the degree of difference reveals the changing pattern of pollution risk, with areas of high degree of difference indicating that the pollution state is undergoing significant changes.

[0082] The spatial location of risk transition zones can be identified based on the degree of difference. Risk transition zones are areas where the pollution risk status is about to change significantly, typically manifested as continuous spatial regions where the degree of difference exceeds a certain threshold. Risk transition zone identification can employ either the region growing method or the contour extraction method. The region growing method starts from the extreme points of difference and gradually incorporates adjacent grid nodes where the degree of difference exceeds the threshold; the contour extraction method directly extracts closed curves where the degree of difference equals the threshold as the boundary of the risk transition zone.

[0083] Obtain the pollution risk level corresponding to the interpolation location, and draw the pollution risk level boundary based on the pollution risk level. The level boundary drawing adopts either the contour line method or the region segmentation method. The contour line method extracts a closed curve where the risk score equals the level boundary threshold; the region segmentation method directly divides the space into regions of different risk levels according to the risk score threshold. Smoothing should be performed during the boundary drawing process to reduce noise and jagged boundaries.

[0084] A spatial distribution map of pollution risk is generated by overlaying the spatial location of the risk transition zone, the boundary of the pollution risk level, and the interpolated pollution risk score. The bottom layer is a continuous color gradient map of the interpolated pollution risk score, the middle layer is the boundary line of the pollution risk level, and the top layer is a highlighted display of the risk transition zone.

[0085] The distribution differences of various pollution risk levels within and outside the risk transition zone are statistically analyzed, including area proportion, spatial clustering, and directionality. Based on these distribution differences, spatial transition patterns of pollution risk levels are identified. Typical patterns include gradient diffusion, transition, and fluctuation. Gradient diffusion is characterized by a continuous and gradual change in risk levels spatially; transition is characterized by a sudden change in risk levels within a specific area; and fluctuation is characterized by the alternating appearance of risk levels in space.

[0086] This process integrates spatial transformation patterns, interpolated pollution evolution trend values, and spatial distribution maps of pollution risks to construct and output a pollution risk assessment report. The report includes an overview of the spatial distribution of pollution risks, analysis of risk transition zones, prediction of risk evolution trends, and management recommendations. Visualizations include multi-scale spatial distribution maps of pollution risks, key analysis maps of risk transition zones, and risk evolution trend maps.

[0087] This invention utilizes a method for generating and outputting pollution risk spatial distribution maps and assessment reports to achieve a spatial visualization and analysis of pollution risks reflected by the environmental microbiome. This method extends pollution risk information from discrete sampling points to continuous space through spatial interpolation technology, overcoming the limitations of traditional point-based assessments and providing a more comprehensive view of regional pollution risks. The identification of risk transition zones and the analysis of spatial transition patterns reveal the spatial dynamics of pollution risks, providing a scientific basis for predicting pollution diffusion trends. The multi-layered information overlay of the risk spatial distribution map visually demonstrates the spatial variation and hierarchical distribution of pollution risks, facilitating environmental managers to quickly identify high-risk and potential-risk areas.

[0088] like Figure 2 As shown, Figure 2 This is a schematic diagram of the structure of a pollution risk assessment system based on environmental microbiome characteristics provided in an embodiment of the present invention. The system includes: The sequencing and annotation module 201 is used to obtain microbial community samples of the target water body at multiple time points for sequencing, and to perform classification annotation and functional annotation on the sequencing results to obtain the abundance of microbial taxonomic units and functional genes. The feature screening module 202 is used to calculate the correlation coefficient between the abundance of microbial taxonomic units and the abundance of functional genes to construct a correlation matrix, and to screen pollution response features based on the structural differences in the correlation matrix between adjacent time points. The risk scoring module 203 is used to obtain the environmental quality standard limit corresponding to the pollution response characteristics, and calculates the weighted deviation of the pollution response characteristics from the environmental quality standard limit through nonlinear mapping as the pollution risk score. The classification module 204 is used to perform time-series fitting and periodic analysis on the pollution risk score to obtain the pollution evolution trend value and pollution cycle parameter. Based on the pollution risk score and pollution evolution trend value, the evolution trajectory is constructed and geometric feature parameters are extracted. The pollution risk level is classified according to the geometric feature parameter and pollution cycle parameter. The spatial analysis module 205 is used to obtain the sampling location coordinates of the microbial community sample as the interpolation location, use the pollution risk score and pollution evolution trend as the interpolation value to perform spatial interpolation, combine the pollution risk level to generate a pollution risk spatial distribution map and output a pollution risk assessment report.

[0089] One technical solution provided in this embodiment of the invention is an electronic device, including: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the steps in any of the aforementioned methods.

[0090] One technical solution provided in this embodiment of the invention is a computer-readable storage medium storing a computer program, wherein the processor executes the computer program to implement the steps in any of the aforementioned methods.

[0091] The specific embodiments described above are preferred embodiments of the present invention and are not intended to limit the specific scope of the present invention. The scope of the present invention includes, but is not limited to, these specific embodiments. All equivalent changes made in accordance with the shape and structure of the present invention are within the protection scope of the present invention.

Claims

1. A pollution risk assessment method based on environmental microbiome characteristics, characterized in that, Includes the following steps: Microbial community samples from the target water body at multiple time points were obtained and sequenced. The sequencing results were then annotated for classification and function to obtain the abundance of microbial taxa and functional genes. The correlation coefficient between the abundance of microbial taxa and the abundance of functional genes was calculated to construct a correlation matrix. Pollution response characteristics were obtained by screening the structural differences in the correlation matrix between adjacent time points. Obtain the environmental quality standard limits corresponding to the pollution response characteristics, and calculate the weighted deviation of the pollution response characteristics from the environmental quality standard limits through nonlinear mapping as the pollution risk score; Time series fitting and periodic analysis are performed on the pollution risk score to obtain the pollution evolution trend value and pollution cycle parameter. Based on the pollution risk score and pollution evolution trend value, the evolution trajectory is constructed and geometric feature parameters are extracted. The pollution risk level is divided according to the geometric feature parameter and pollution cycle parameter. The sampling location coordinates of the microbial community samples are obtained as the interpolation location. Spatial interpolation is performed using the pollution risk score and pollution evolution trend value as the interpolation values. Combined with the pollution risk level, a spatial distribution map of pollution risk is generated and a pollution risk assessment report is output.

2. The method according to claim 1, characterized in that, Microbial community samples from the target water body at multiple time points were obtained and sequenced. The sequencing results were then annotated for taxonomic and functional data, yielding the abundance of microbial taxa and functional genes, including: Water samples were collected from the target water body at multiple time points as microbial community samples, and nucleic acids were extracted from the microbial community samples for high-throughput sequencing to obtain sequencing results. Sequence clustering is performed on the sequencing results. Sequences with similarity exceeding a preset similarity threshold are grouped into the same microbial taxonomic unit to complete the classification annotation. The abundance of the initial microbial taxonomic unit is obtained by calculating the ratio of the number of sequences in each microbial taxonomic unit to the total number of sequences in the sequencing results. The frequency of occurrence of each microbial taxonomic unit at multiple time points was statistically analyzed, and the initial microbial taxonomic unit abundance corresponding to the microbial taxonomic units with an occurrence frequency higher than a preset frequency threshold was selected and normalized to obtain the microbial taxonomic unit abundance. The functional gene types encoded by each sequence in the sequencing results are identified. The metabolic pathways to which the functional gene types belong are grouped. The abundance of the metabolic pathway is obtained by summing the proportion of the number of sequences contained in each metabolic pathway to the total number of sequences. The functional gene abundance is obtained by normalization.

3. The method according to claim 1, characterized in that, A correlation matrix was constructed by calculating the correlation coefficient between the abundance of microbial taxa and the abundance of functional genes. Pollution response characteristics were then selected based on the structural differences in the correlation matrix between adjacent time points. For each time point, the time series correlation coefficient between the abundance of each microbial taxonomic unit and the abundance of each functional gene is calculated to obtain the correlation coefficient. All correlation coefficients are then filled into a correlation matrix according to the correspondence between microbial taxonomic units and functional genes. Singular value decomposition is performed on the correlation matrix at each time point to obtain the dominant singular values ​​and the corresponding dominant singular vectors. Vector projection is performed on the dominant singular vectors of adjacent time points to calculate the projection residuals. The projection residuals and the changes in the dominant singular values ​​are jointly weighted to obtain the structural difference values ​​of the correlation matrix between adjacent time points, and a sequence of structural difference values ​​is constructed. Baseline fitting is performed on the structural difference value sequence to identify structural anomaly time points that deviate from the baseline. The correlation matrix corresponding to the structural anomaly time point is extracted and the correlation matrix of the baseline time period is calculated to obtain the difference matrix. The difference matrix is ​​divided into blocks and clustered according to the numerical distribution of matrix elements to obtain response modules. The concentration of matrix elements within each response module is calculated, and the microbial taxonomic unit and functional gene corresponding to the response module with the highest concentration are extracted as pollution response features.

4. The method according to claim 1, characterized in that, The environmental quality standard limits corresponding to pollution response characteristics are obtained, and the weighted deviation of the pollution response characteristics from the environmental quality standard limits is calculated using a nonlinear mapping to obtain the pollution risk score, which includes: Functional genes are extracted from pollution response characteristics, the metabolic pathways to which the functional genes belong are identified, the integrity of the metabolic pathways is calculated, pollutant types are screened based on the integrity of the metabolic pathways, and the environmental quality standard limits corresponding to the pollutant types are obtained. The abundance of microbial taxa and functional genes corresponding to pollution response characteristics are extracted. The normalized ratios of microbial taxa abundance and functional gene abundance relative to environmental quality standard limits are calculated respectively. The deviation of pollution response characteristics from environmental quality standard limits is obtained by nonlinear mapping transformation of the normalized ratios. Based on the correlation coefficients of pollution response characteristics in the correlation matrix, the contribution rate of the correlation coefficients is calculated to construct a weight vector. The deviation is then weighted and summed using the weight vector to obtain the weighted deviation as the pollution risk score.

5. The method according to claim 1, characterized in that, Time-series fitting and periodic analysis are performed on pollution risk scores to obtain pollution evolution trend values ​​and pollution cycle parameters. Based on the pollution risk scores and pollution evolution trend values, an evolution trajectory is constructed and geometric feature parameters are extracted. Pollution risk levels are then classified according to the geometric feature parameters and pollution cycle parameters, including: A time series of pollution risk scores is constructed, the time series is segmented by a sliding window, the average rate of change within each segment is calculated, abrupt changes in the average rate of change between adjacent segments are identified, and the time period is re-divided based on the abrupt changes. The weighted combination of the average rate of change in each time period is calculated as the pollution evolution trend value. Wavelet decomposition is performed on the time series to extract the periodic components, and the dominant oscillation period of the periodic components is calculated as the pollution period parameter. The pollution risk scores are arranged in chronological order as a state vector, and the pollution evolution trend values ​​are arranged in chronological order as a driving vector. The projection trajectory of the state vector in the direction of the driving vector is calculated, and the oscillation amplitude and phase drift of the projection trajectory are extracted as geometric feature parameters through morphological analysis. Risk fluctuation intensity is constructed based on the oscillation amplitude in the geometric characteristic parameters, and periodic stability index is constructed based on the pollution cycle parameters. The coupling index is calculated based on the risk fluctuation intensity and periodic stability index. Multi-level thresholds are set according to the numerical distribution of the coupling index, and the samples are divided into different pollution risk levels according to the multi-level thresholds.

6. The method according to claim 5, characterized in that, Calculate the projected trajectory of the state vector in the direction of the driving vector, and perform morphological analysis on the projected trajectory to extract the oscillation amplitude and phase drift of the trajectory, including: Calculate the unit direction vector of the driving vector, calculate the dot product of the state vector and the unit direction vector to obtain the projection scalar sequence, and connect the projection scalar sequence in time order to form the projection trajectory. The upper and lower envelopes are obtained by fitting the projection trajectory with an envelope. The instantaneous distance sequence between the upper and lower envelopes is calculated. The oscillation amplitude is obtained by performing a time-weighted average on the instantaneous distance sequence. The extreme points of the projected trajectory are detected to obtain peak points and valley points. The time interval between adjacent peak points is calculated to obtain the peak interval sequence. The average value of the peak interval sequence is calculated, and the difference between the average value and the pollution cycle parameter is used as the phase drift.

7. The method according to claim 1, characterized in that, The sampling location coordinates of the microbial community samples are obtained as the interpolation location. Spatial interpolation is performed using the pollution risk score and pollution evolution trend value as interpolation values. Combined with the pollution risk level, a spatial distribution map of pollution risk is generated, and a pollution risk assessment report is output, including: The sampling location coordinates of the microbial community sample are obtained as the interpolation location, and the interpolation location is divided into spatial grids to obtain grid nodes; Calculate the spatial distance between grid nodes and interpolation locations, construct distance weights based on the spatial distances, use pollution risk scores and pollution evolution trend values ​​as interpolation values, and use the distance weights to perform spatial interpolation on the interpolation values ​​to obtain the interpolated pollution risk scores and interpolated pollution evolution trend values ​​for each grid node. Calculate the difference between the interpolated pollution risk score and the interpolated pollution evolution trend value for each grid node, and identify the spatial location of the risk transition zone based on the difference. Obtain the pollution risk level corresponding to the interpolation location, draw the pollution risk level boundary according to the pollution risk level, and generate a pollution risk spatial distribution map by superimposing the spatial location of the risk transition zone, the pollution risk level boundary and the interpolated pollution risk score. The distribution differences of each pollution risk level inside and outside the risk transition zone are statistically analyzed. Based on the distribution differences, the spatial transition patterns of pollution risk levels are identified. The spatial transition patterns, interpolated pollution evolution trend values, and pollution risk spatial distribution maps are integrated to construct and output a pollution risk assessment report.

8. A pollution risk assessment system based on environmental microbiome characteristics, used to implement the method described in any one of claims 1-7, characterized in that, The system includes: The sequencing and annotation module is used to obtain microbial community samples of the target water body at multiple time points for sequencing, and to perform classification annotation and functional annotation on the sequencing results to obtain the abundance of microbial taxonomic units and functional genes. The feature screening module is used to calculate the correlation coefficient between the abundance of microbial taxonomic units and the abundance of functional genes to construct a correlation matrix, and to screen pollution response features based on the structural differences in the correlation matrix between adjacent time points. The risk scoring module is used to obtain the environmental quality standard limits corresponding to the pollution response characteristics. It calculates the weighted deviation of the pollution response characteristics from the environmental quality standard limits through nonlinear mapping as the pollution risk score. The classification module is used to perform time-series fitting and periodic analysis on pollution risk scores to obtain pollution evolution trend values ​​and pollution cycle parameters. Based on the pollution risk scores and pollution evolution trend values, the module constructs evolution trajectories and extracts geometric feature parameters. Based on the geometric feature parameters and pollution cycle parameters, the module classifies pollution risk levels. The spatial analysis module is used to obtain the sampling location coordinates of microbial community samples as interpolation locations, use pollution risk scores and pollution evolution trend values ​​as interpolation values ​​for spatial interpolation, combine pollution risk levels to generate a pollution risk spatial distribution map and output a pollution risk assessment report.

9. An electronic device, characterized in that, include: A memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor, when executing the computer program, implements the steps of the method as described in any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores computer program instructions that, when executed by a processor, implement the steps of the method as described in any one of claims 1 to 7.

Images