Method for quantifying influence of climate change and human activity on watershed vegetation

Abstract

The invention discloses a method for quantifying the influence of climate change and human activity on watershed vegetation, which comprises the following steps: acquiring data of a research area from a plurality of data sources, performing time-space consistency processing, and calculating saturation vapor pressure difference; performing feature extraction on the processed data, and constructing a high-dimensional feature tensor; constructing a space-time diagram based on the high-dimensional feature tensor, performing dimension reduction and dynamic modeling by applying a map embedding technology, and analyzing space-time dynamic features of watershed vegetation; separating and quantifying direct influence of climate and human activities on vegetation by adopting a causal network and a regression decomposition method; comprehensive evaluation is carried out, and uncertainty analysis is carried out. According to the method, the ability of understanding the dynamic state of a complex ecosystem is remarkably improved, powerful scientific support is provided for watershed ecosystem management, and more accurate, reliable and comprehensive decision support can be provided.

Description

Technical Field

[0001] The invention relates to the field of quantification technology, in particular to a method for quantifying the impact of climate change and human activities on watershed vegetation. Background Art

[0002] Against the backdrop of global climate change and increasing human activities, accurately quantifying their impacts on watershed vegetation has become a key research topic in ecology and environmental science. As an important ecological unit, the vegetation coverage of a watershed directly affects soil and water conservation, carbon cycle, and biodiversity. Therefore, a deep understanding of the impact mechanisms of climate change and human activities on watershed vegetation will not only help assess the vulnerability and resilience of ecosystems, but also provide a scientific basis for formulating effective ecological management strategies and climate change adaptation measures. Especially in transboundary watersheds, such research has important practical significance for coordinating upstream and downstream relations and promoting regional sustainable development.

[0003] At present, researchers have developed a variety of methods to assess the impact of climate change and human activities on vegetation. Common methods include vegetation index analysis based on remote sensing data, statistical modeling, and ecosystem model simulation. In terms of data acquisition, the integration of multi-source remote sensing data makes large-scale vegetation monitoring possible. In terms of data analysis, researchers usually use statistical methods such as correlation analysis, regression analysis, and trend analysis to explore the relationship between vegetation changes and environmental factors. Some studies have also tried to use machine learning algorithms, such as random forests and support vector machines, to improve the predictive ability of the model. In terms of impact separation, residual analysis is widely used to distinguish the relative contributions of climate change and human activities. In addition, some researchers also use ecosystem models, such as BIOME-BGC or LPJ-GUESS, to simulate the response of vegetation to environmental changes.

[0004] However, the current research methods still have some limitations and technical challenges. First, in the data preprocessing stage, the spatiotemporal consistency processing of multi-source data is still a thorny issue, especially when dealing with data of different spatial resolutions and temporal frequencies. Second, traditional linear models have difficulty capturing the nonlinear response of vegetation to environmental factors, especially in the presence of ecological thresholds. Furthermore, existing impact separation methods often ignore the interactions between factors, which may lead to incorrect estimates of the impact of a single factor. In terms of time scale, most studies focus on changes on interannual or longer time scales, ignoring intra-seasonal and inter-annual variability. Spatially, how to effectively deal with the spatial heterogeneity and scale dependence of vegetation change remains a challenge. In addition, in causal analysis, traditional correlation analysis has difficulty distinguishing true causality from spurious correlation. Finally, existing uncertainty analysis methods are often oversimplified and it is difficult to fully evaluate the reliability of prediction results, especially in the face of complex ecosystems and multiple future scenarios. Solving these technical problems is crucial to improving the accuracy and reliability of studies on watershed vegetation changes. Summary of the invention

[0005] In order to solve the problems existing in the above-mentioned background technology, the present invention is proposed.

[0006] Specifically, the present invention provides the following technical solutions: An embodiment of the present invention provides a method for quantifying the impact of climate change and human activities on watershed vegetation, which includes obtaining data of a study area from multiple data sources, performing spatiotemporal consistency processing to obtain a preprocessed data set in a unified format, and calculating the saturated water vapor pressure difference; extracting features based on the preprocessed data set, and fusing the features with the preprocessed data set to construct a high-dimensional feature tensor; constructing a spatiotemporal graph based on the high-dimensional feature tensor, applying graph embedding technology to reduce the dimension and establish a dynamic response model, and analyzing the spatiotemporal dynamic characteristics of watershed vegetation; using causal networks and regression decomposition methods to separate and quantify the direct impact of climate and human activities on vegetation; and conducting a comprehensive evaluation and uncertainty analysis.

[0007] The beneficial effects of the present invention are as follows: the present invention significantly improves the ability to understand the dynamics of complex ecosystems, and provides strong scientific support for watershed ecosystem management. In particular, when facing the challenges of watershed ecosystem management under the background of global change, the present invention can provide more accurate, reliable and comprehensive decision-making support, which can not only assess the current ecological status, but also predict future trends, and provide a scientific basis for the formulation of long-term ecological protection and restoration strategies. In addition, the flexibility and scalability of the present invention make it suitable for the study of watershed ecosystems of different types and scales, and provide new tools and ideas for comparative studies between regions and ecological assessments at the global scale. BRIEF DESCRIPTION OF THE DRAWINGS

[0008] Figure 1 It is a flow chart of the present invention.

[0009] Figure 2 It is a flowchart for extracting features, fusing the features with processed data, and constructing a high-dimensional feature tensor.

[0010] Figure 3 It is a flow chart of applying the improved Markov chain model to extract the characteristics of land use change.

[0011] Figure 4 It is a flow chart that uses nonlinear response functions to quantify the effect of carbon dioxide fertilization based on carbon dioxide concentration data and vegetation index data.

[0012] Figure 5 It is a flowchart for constructing a high-dimensional feature tensor. DETAILED DESCRIPTION

[0013] For a detailed description of the implementation process of this program, refer to Figures 1 to 5 , the specific description is as follows: According to one aspect of the present application, the calculation of the saturated water vapor pressure difference includes: Calculation of saturated water vapor pressure using the modified Magnus formula , the formula is: ; Wherein, T is the temperature data; Using the calculated saturated water vapor pressure and relative humidity data , calculate the actual water vapor pressure, the formula is: ; Among them, RH is the relative humidity data; Calculating Saturation Vapor Pressure Deficit , the formula is: VPD = es - ea.

[0014] According to one aspect of the present application, the extraction of features and fusion of the features with processed data to construct a high-dimensional feature tensor includes the following steps: The MTPI algorithm is used to extract terrain features based on multi-scale windows; Applying the improved Markov chain model to extract the characteristics of land use change; Based on CO2 concentration data and vegetation index data, a nonlinear response function is used to quantify the CO2 fertilization effect. The processed data are fused with the extracted terrain features, land use change features, and CO2 fertilization effect features to obtain a high-dimensional feature tensor.

[0015] According to one aspect of the present application, the application of the improved Markov chain model to extract land use change characteristics comprises the following steps: Extracting land use change data from the preprocessed data set and storing it as a land use time series data set; For every two adjacent time points in the land use time series data set, the transition frequency between the land use types is calculated, the transition frequency is divided by the total frequency of the initial state to obtain a transition probability matrix P, and the calculated transition probability is stored as a transition probability matrix data set; For each spatial pixel, using the probability in the transition probability matrix data set, repeatedly simulate the land use change sequence within the study period to generate multiple change sequences, and store the generated change sequences as a land use change trajectory data set; Applying a sequence clustering algorithm to the land use change trajectory dataset, determining an optimal number of clusters, assigning a cluster label to each spatial pixel, and storing the clustering results as a land use change pattern dataset; Merging the land use change pattern dataset with the land use time series dataset, calculating the statistical characteristics of each change pattern, and storing the merged data as a land use change characteristic dataset; The land use time series dataset is a collection of extracted land use change data.

[0016] According to one aspect of the present application, the method of quantifying the carbon dioxide fertilization effect using a nonlinear response function based on carbon dioxide concentration data and vegetation index data comprises the following steps: Extracting carbon dioxide concentration data and normalized vegetation index data from the preprocessed data set, and storing them as a carbon dioxide concentration data set and an NDVI data set, respectively; Construct a nonlinear response function model, where the nonlinear response function is: , Among them, a and b are the parameters to be estimated. is the carbon dioxide concentration, is the reference value of carbon dioxide concentration; For each spatial pixel, a nonlinear least square method is used, with the logarithmic concentration ratio data set as an independent variable and the NDVI data set as a dependent variable, to fit the nonlinear response function, and the estimated a and b parameters are stored as a parameter estimation data set; Using the parameters in the parameter estimation dataset, calculate the predicted NDVI value for each spatiotemporal location; The parameters in the parameter estimation dataset and the logarithmic concentration ratio dataset are used to calculate the carbon dioxide fertilization effect at each spatiotemporal location, and the calculation results are stored as a carbon dioxide fertilization effect dataset.

[0017] According to one aspect of the present application, the construction process of the high-dimensional feature tensor is as follows: Aligning the preprocessing dataset, the terrain feature dataset, the land use change feature dataset, and the carbon dioxide fertilization effect dataset in time and space dimensions, and storing them as an aligned feature dataset; reorganizing the aligned feature dataset into a three-dimensional tensor, and storing the constructed tensor as an original feature tensor dataset; Performing Tucker decomposition on the original feature tensor data set, initializing the core tensor and factor matrix, and iteratively optimizing the decomposition result using alternating least squares method until convergence or reaching a maximum number of iterations; Use the cross-validation method to determine the optimal decomposition rank and store it as a decomposition rank data set; Using the optimal rank in the decomposition rank dataset, performing a final Tucker decomposition on the original feature tensor dataset, and storing the obtained core tensor and factor matrix as a decomposition result dataset; The high-dimensional feature tensor is reconstructed using the core tensor and factor matrix in the decomposition result data set, and the reconstructed tensor is flattened into a two-dimensional matrix, in which each row represents a spatiotemporal position and each column represents a fused feature, and the two-dimensional matrix is ​​stored as a high-dimensional feature data set.

[0018] According to one aspect of the present application, constructing a spatiotemporal graph based on the high-dimensional feature tensor includes: Extract spatial coordinates and time information from the high-dimensional feature data set, and construct a spatiotemporal weighted graph G(V,E), where V is a vertex set and E is an edge set. Take each spatiotemporal position as a vertex of the spatiotemporal weighted graph and calculate the weight of the edge: , in, is the spatial distance, and is the time point, σ and τ are adjustable parameters; The optimal σ and τ values ​​are determined by cross-validation, and a weighted adjacency matrix is ​​constructed based on the calculated weights to obtain a spatiotemporal weighted graph data structure.

[0019] According to one aspect of the present application, the causal network and regression decomposition method are used to separate and quantify the gas Direct impacts of climate and human activities on vegetation include: Use causal entropy method to construct causal network and quantify the direct causal relationship between climate and human activities; Separate and quantify the contribution of climate and human activities to vegetation through regression decomposition methods; The generalized additive model was used to quantify the nonlinear effects and thresholds of each influencing factor; The dynamic time warping algorithm is used to analyze the time lag effect and generate the lag effect matrix.

[0020] Embodiment 1: This embodiment provides a method for quantifying the impact of climate change and human activities on watershed vegetation, comprising the following steps: Step S1, obtaining vegetation remote sensing data, climate data and auxiliary data of the study area from multiple data sources, performing spatiotemporal consistency processing on the obtained data, and calculating the saturated water vapor pressure difference.

[0021] First, the normalized vegetation index time series data, vegetation functional type data and vegetation coverage data of the study area were obtained from the remote sensing database. Then, the precipitation data, temperature data and relative humidity data of the study area were obtained from the meteorological database. Next, the digital elevation model data, land use change data and carbon dioxide concentration data of the study area were obtained from the geographic information database. All the acquired data were interpolated in time series, resampled in spatial resolution and projected to obtain a preprocessed data set in a unified format. Finally, the saturated water vapor pressure deficit data were calculated using the improved Magnus formula using the preprocessed temperature data and relative humidity data.

[0022] Specifically, step S11 includes the following steps: Step S11, obtain the normalized vegetation index time series data of the study area from the remote sensing database. Then, obtain the vegetation function type data of the study area from the same database. Then, obtain the vegetation coverage data of the study area from the database. The three types of data obtained are stored as the normalized vegetation index data set, the vegetation function type data set and the vegetation coverage data set respectively.

[0023] Step S12: Obtain precipitation data of the study area from the meteorological database. Next, obtain temperature data of the study area from the same database. Then, obtain relative humidity data of the study area from the database. The three types of data obtained are stored as precipitation data set, temperature data set and relative humidity data set respectively.

[0024] Step S13, obtaining digital elevation model data of the study area from the geographic information database. Next, obtaining land use change data of the study area from the same database. Then, obtaining carbon dioxide concentration data of the study area from the database. The three types of data obtained are stored as digital elevation model datasets, land use change datasets, and carbon dioxide concentration datasets, respectively.

[0025] Step S14, perform time series interpolation processing on the normalized vegetation index dataset, vegetation functional type dataset and vegetation coverage dataset obtained in step S11. Then, perform time series interpolation processing on the precipitation dataset, temperature dataset and relative humidity dataset obtained in step S12. Next, perform time series interpolation processing on the digital elevation model dataset, land use change dataset and carbon dioxide concentration dataset obtained in step S13. Resample the spatial resolution of all interpolated datasets to have the same spatial resolution. Then, perform projection transformation on all resampled datasets to have the same map projection. Finally, merge all processed datasets to obtain a preprocessed dataset in a unified format.

[0026] Step S15, extracting temperature data and relative humidity data from the preprocessed data set obtained in step S14. Calculate the saturated water vapor pressure using the improved Magnus formula, the formula is: es = 0.611 * exp((17.27 * T) / (T +237.3)), where T is the temperature data. Then, use the calculated saturated water vapor pressure and relative humidity data to calculate the actual water vapor pressure, the formula is: ea = es * (RH / 100), where RH is the relative humidity data. Finally, calculate the saturated water vapor pressure difference, the formula is: VPD = es - ea. Store the calculated saturated water vapor pressure difference data as a saturated water vapor pressure difference data set.

[0027] Through multi-source data acquisition and preprocessing, comprehensive data coverage of climate change and human activities in the study area is achieved. This method not only integrates remote sensing data, climate data and geographic information data, but also introduces the calculation of saturated vapor pressure deficit (VPD), which is an innovation of traditional methods. In particular, in the preprocessing stage, time series interpolation, spatial resolution resampling and projection transformation techniques are used to ensure the spatiotemporal consistency of data from different sources. This preprocessing method solves the problem of mismatch between multi-source data in time and space scales, and provides high-quality input data for subsequent analysis. The introduction of VPD calculation is an important supplement to the traditional method that only considers temperature and precipitation, because VPD directly affects the transpiration and water use efficiency of plants, and can more accurately reflect the impact of climate change on vegetation. This comprehensive and sophisticated data processing method provides a solid data foundation for quantifying the impact of climate change and human activities on vegetation in the basin, and significantly improves the accuracy and reliability of subsequent analysis.

[0028] Step S2, using the data obtained and preprocessed in step S1, extract terrain features, land use change features and carbon dioxide fertilization effect features, and fuse these features with the preprocessed data. First, based on the digital elevation model data in the preprocessed data set, a multi-scale terrain position index algorithm is used to extract terrain features. Then, an improved Markov chain model is applied to the land use change data in the preprocessed data set to extract land use change features. Next, based on the carbon dioxide concentration data and normalized vegetation index data in the preprocessed data set, a nonlinear response function is used to quantify the carbon dioxide fertilization effect. Finally, the preprocessed data set is fused with the extracted terrain features, land use change features and carbon dioxide fertilization effect features using a tensor decomposition method to obtain a high-dimensional feature tensor.

[0029] Specifically, step S2 is further as follows: Step S21, extracting digital elevation model data from the preprocessed data set obtained in step S14. Applying a multi-scale terrain position index algorithm to the extracted digital elevation model data. First, set multiple spatial scales, such as 3x3, 5x5, and 7x7 pixel windows. For each scale, calculate the elevation difference between the central pixel and the surrounding pixels. Then, perform a weighted average on the elevation difference of each scale, with the weight set inversely proportional to the distance. Next, normalize the weighted average results of different scales. Finally, merge the normalized multi-scale terrain position index results to obtain a terrain feature data set.

[0030] Step S22, extracting land use change data from the preprocessed data set obtained in step S14. Applying the improved Markov chain model to the extracted land use change data. First, calculating the transition probability matrix between each land use type between adjacent time steps. Then, based on the transition probability matrix, calculating the land use change trajectory of each pixel during the study period. Next, performing cluster analysis on the calculated change trajectory to identify the main land use change patterns. Finally, combining the clustering results with the original land use change data to obtain a land use change feature dataset.

[0031] Specifically, step S22 is further as follows: Step S221: extract land use change data from the preprocessed data set obtained in step S14, and store it as a land use time series data set.

[0032] Step S222: Construct a land use type transition matrix. For every two adjacent time points in the land use time series data set, calculate the transition frequency between each land use type. Divide the transition frequency by the total frequency of the initial state to obtain the transition probability. Store the calculated transition probability as a transition probability matrix data set.

[0033] Step S223, calculate the land use change trajectory. For each spatial pixel, use the probability in the transition probability matrix data set to simulate its land use change sequence during the study period. Repeat this process multiple times to generate multiple possible change sequences. Store the generated change sequence as a land use change trajectory data set.

[0034] Step S224: Perform cluster analysis of the change trajectory. Apply a sequence clustering algorithm, such as dynamic time warping clustering, to the land use change trajectory dataset. Determine the optimal number of clusters and assign a cluster label to each spatial pixel. Store the clustering results as a land use change pattern dataset.

[0035] Step S225: Combine the original land use change data with the clustering results. Merge the land use change pattern dataset with the land use time series dataset in step S221. Calculate the statistical characteristics of each change pattern, such as duration, change frequency, etc. Store the merged data as a land use change feature dataset.

[0036] First, the land use type transition matrix is ​​constructed, which not only captures the direction of land use change, but also quantifies the probability of change, which provides important information for understanding the dynamic process of land use change. Secondly, by calculating the land use change trajectory, this method can describe the historical change sequence of each spatial unit. This time series analysis method can reveal the long-term land use change pattern. Thirdly, the sequence clustering algorithm is applied to perform change trajectory clustering analysis. This method can identify typical land use change patterns and help understand the land use change trend at the regional scale. Finally, by combining the original land use change data with the clustering results, the statistical characteristics of each change pattern, such as duration and change frequency, are calculated, which provides detailed indicators for quantifying the impact of human activities on vegetation. In the study of watershed vegetation, this improved land use change analysis method can accurately characterize the impact process and intensity of human activities on vegetation. It can not only identify key land use change types (such as urbanization, agricultural expansion or deforestation), but also quantify the spatiotemporal patterns and intensity of these changes. This detailed information is crucial to understanding how human activities change vegetation cover, affect ecosystem functions, and interact with climate change. In addition, the output of this method can be directly used to parameterize ecological models and improve the accuracy of the models' simulations of the impacts of human activities.

[0037] Step S23, extracting carbon dioxide concentration data and normalized vegetation index data from the preprocessed data set obtained in step S14. Using a nonlinear response function to quantify the carbon dioxide fertilization effect. First, set a reference carbon dioxide concentration value. Then, for each pixel, calculate the natural logarithm of the ratio of the actual carbon dioxide concentration to the reference concentration. Then, substitute the calculated result into the nonlinear response function , where a and b are fitting parameters. The optimal a and b parameters for each pixel are determined by minimizing the error between the actual NDVI and the model prediction value. Finally, the CO2 fertilization effect is calculated using the determined parameters to obtain the CO2 fertilization effect dataset.

[0038] Specifically, step S23 is further as follows: Step S231, extracting carbon dioxide concentration data and normalized difference vegetation index data from the preprocessed data set obtained in step S14, and storing them as a carbon dioxide concentration data set and a NDVI data set respectively.

[0039] Step S232: Set a reference carbon dioxide concentration value. Calculate the natural logarithm of the ratio of the carbon dioxide concentration at each spatiotemporal position to the reference concentration, and store it as a logarithmic concentration ratio data set.

[0040] Step S233: construct a nonlinear response function model. Define the function , where a and b are the parameters to be estimated.

[0041] Step S234: For each spatial pixel, use the nonlinear least squares method, with the logarithmic concentration ratio data set as the independent variable and the NDVI data set as the dependent variable, to fit the nonlinear response function in S233. The estimated a and b parameters are stored as a parameter estimation data set.

[0042] Step S235: Calculate the predicted NDVI value for each spatiotemporal location using the parameters in the parameter estimation dataset. Calculate the root mean square error and the coefficient of determination between the predicted value and the actual value in the NDVI dataset. Store the evaluation results as a model evaluation dataset.

[0043] Step S236: Calculate the carbon dioxide fertilization effect at each spatiotemporal position using the parameters in the parameter estimation data set and the logarithmic concentration ratio data set, and store the calculation results as a carbon dioxide fertilization effect data set.

[0044] First, the nonlinear response function model can more accurately describe the effect of vegetation on the atmosphere. The complex response of concentration changes is closer to the actual ecological process than the traditional linear model. concentration and calculate the logarithmic concentration ratio, which takes into account the effect of vegetation on The sensitivity of relative concentration changes, rather than absolute concentration, is more consistent with plant physiological principles. Thirdly, the nonlinear least squares method is used for parameter estimation, and the parameters are fitted separately for each spatial pixel. This method takes into account spatial heterogeneity and can capture the sensitivity of vegetation in different regions to changes in concentration. Finally, through the model evaluation step, the error and determination coefficient between the predicted value and the actual value are calculated to ensure the reliability of the quantitative results. The quantification method of fertilization effect is of great significance. It can separate the The changes in vegetation productivity caused by increased concentrations are crucial for accurately assessing the impact of climate change on vegetation. Especially in long-term ecosystem studies, this method can help researchers distinguish The fertilization effect is combined with the influence of other environmental factors such as temperature and precipitation changes. In addition, because the method takes into account spatial heterogeneity, it is able to reveal the impact of different types of vegetation or different growing conditions on The differences in sensitivity to changes in concentration provide an important basis for predicting vegetation responses under future climate change. Quantification of fertilization effects has important application value in improving global carbon cycle models, assessing ecosystem carbon sequestration potential, and formulating relevant climate change adaptation strategies.

[0045] Step S24, fuse the preprocessed data set obtained in step S14 with the terrain feature data set obtained in step S21, the land use change feature data set obtained in step S22, and the carbon dioxide fertilization effect data set obtained in step S23. First, align all data sets in time and space dimensions. Then, construct a high-dimensional data tensor, in which each dimension corresponds to space, time, and features. Next, apply a tensor decomposition method, such as Tucker decomposition or CP decomposition, to the constructed high-dimensional data tensor. During the decomposition process, the optimal decomposition rank is determined by cross-validation. Then, reconstruct the data tensor based on the decomposition results to obtain a fused high-dimensional feature tensor. Finally, flatten the high-dimensional feature tensor into a two-dimensional matrix, in which each row represents a spatiotemporal position and each column represents a fused feature, to obtain a high-dimensional feature data set.

[0046] Specifically, step S24 is further as follows: Step S241, align the preprocessed data set obtained in step S14, the terrain feature data set obtained in step S21, the land use change feature data set obtained in step S225, and the carbon dioxide fertilization effect data set obtained in step S236 in time and space dimensions, and store them as an aligned feature data set.

[0047] Step S242: reorganize the aligned feature data set into a three-dimensional tensor, wherein the first dimension represents space, the second dimension represents time, and the third dimension represents features. The constructed tensor is stored as an original feature tensor data set.

[0048] Step S243: Perform Tucker decomposition on the original feature tensor data set, initialize the core tensor and factor matrix, and iteratively optimize the decomposition result using the alternating least squares method until convergence or the maximum number of iterations is reached.

[0049] Step S244: Use the cross-validation method to try different decomposition rank combinations. For each combination, calculate the reconstruction error and the explained variance ratio. Select the optimal rank that can strike a balance between the reconstruction error and the model complexity. Store the optimal rank as a decomposition rank data set.

[0050] Step S245: Use the optimal rank in the decomposition rank dataset to perform a final Tucker decomposition on the original feature tensor dataset, and store the obtained core tensor and factor matrix as a decomposition result dataset.

[0051] Step S246: Use the core tensor and factor matrix in the decomposition result data set to reconstruct the high-dimensional feature tensor. Flatten the reconstructed tensor into a two-dimensional matrix, where each row represents a spatiotemporal position and each column represents a fused feature. Store this matrix as a high-dimensional feature data set.

[0052] In this embodiment, first, multiple environmental features (such as topography, land use change, Fertilization effect, etc.) are organized into high-dimensional data tensors, retaining the original spatiotemporal structure of the data. This method can better capture the complex relationship between data than traditional matrix analysis. Secondly, advanced tensor decomposition techniques such as Tucker decomposition or CP decomposition are used. These methods can effectively process high-dimensional data and extract the most important feature combinations. Thirdly, the optimal decomposition rank is determined by cross-validation. This method can find the best balance between data compression and information retention, avoiding the problem of overfitting or underfitting. Finally, the high-dimensional feature tensor is reconstructed into a two-dimensional matrix, so that each spatiotemporal position has a comprehensive feature representation, which provides high-quality input data for subsequent analysis and modeling. In the study of watershed vegetation, this feature fusion method has important application value. It can effectively integrate environmental data from different sources and scales, overcoming the problems of high data dimension, information redundancy, and strong correlation between features in traditional methods. Through dimensionality reduction and feature extraction, this method greatly reduces the complexity of the data while retaining key ecological information. This not only improves the efficiency and accuracy of subsequent modeling, but also reveals the potential interactions between different environmental factors. Especially when studying the combined impacts of climate change and human activities on vegetation, this approach can simultaneously consider the combined effects of multiple factors, providing a more comprehensive and systematic analytical framework.

[0053] Step S3, based on the high-dimensional feature data set obtained in step S2, construct a spatiotemporal weighted graph, then perform graph embedding, and finally establish a vegetation dynamic response model. First, a spatiotemporal weighted graph is constructed using the high-dimensional feature data set, where vertices represent spatial positions and edges represent spatiotemporal relationships. Then, a spectral clustering algorithm is applied to the constructed graph to obtain a low-dimensional spatiotemporal feature representation. Finally, based on the low-dimensional spatiotemporal feature representation, a nonlinear autoregressive exogenous model is used to construct a vegetation dynamic response model.

[0054] Through advanced feature extraction and fusion, comprehensive consideration and in-depth analysis of complex environmental factors are achieved. The innovation of this step is reflected in several aspects: First, the multi-scale terrain position index (MTPI) algorithm is used to extract terrain features, which can more comprehensively characterize the impact of terrain on vegetation distribution compared to the traditional single-scale method. Secondly, the use of an improved Markov chain model to analyze land use changes can not only capture the trend of changes, but also identify the pattern and frequency of changes, which is crucial to understanding the long-term impact of human activities on vegetation. Thirdly, the introduction of nonlinear response functions to quantify Fertilization effect, which is an important consideration for vegetation response in the context of global change. Finally, a tensor decomposition method is used for feature fusion, which can effectively process high-dimensional data and retain the spatiotemporal structural information of the original data. This multi-angle and multi-scale feature extraction and fusion method enables the analysis to simultaneously consider the impact of topography, land use change, and atmospheric composition change on vegetation, greatly improving the ability to understand the dynamics of complex ecosystems. Especially at the watershed scale, this method can capture the interaction between different environmental factors and provide a powerful tool for comprehensively assessing the combined impact of climate change and human activities on vegetation.

[0055] Step S31, extracting spatial coordinates and time information from the high-dimensional feature data set obtained in step S24. Construct a spatiotemporal weighted graph based on the extracted information. First, take each spatiotemporal position as a vertex of the graph. Then, calculate the spatial distance and time difference between each pair of vertices. Next, use the formula w_ij = exp(-d_ij^2 / σ^2) * exp(-|t_i -t_j| / τ) to calculate the weight of the edge, where d_ij is the spatial distance, t_i and t_j are time points, and σ and τ are adjustable parameters. Determine the optimal σ and τ values ​​by cross-validation. Finally, construct a weighted adjacency matrix based on the calculated weights to obtain a spatiotemporal weighted graph data structure.

[0056] Step S32, apply the spectral clustering algorithm to the spatiotemporal weighted graph constructed in step S31. First, calculate the degree matrix of the graph based on the weighted adjacency matrix. Then, calculate the Laplace matrix L = D - W of the graph, where D is the degree matrix and W is the weighted adjacency matrix. Next, solve the generalized eigenvalue problem L * φ_k = λ_k * D * φ_k to obtain the eigenvalues ​​and eigenvectors. Then, select the eigenvectors corresponding to the K smallest non-zero eigenvalues ​​to form a low-dimensional embedding matrix. Finally, take the coordinates of each vertex in the low-dimensional embedding space as its new feature representation to obtain a low-dimensional spatiotemporal feature data set.

[0057] Specifically, step S32 is further: Step S321, extracting a weighted adjacency matrix from the spatiotemporal weighted graph data structure obtained in step S31, and storing it as a weighted adjacency matrix data set.

[0058] Step S322: Calculate the degree matrix of the graph. Sum each row in the weighted adjacency matrix data set to obtain the degree of each vertex. Store the calculation result as a degree matrix data set.

[0059] Step S323: Calculate the Laplacian matrix of the graph. Use the formula L = D - W, where D is the degree matrix in the degree matrix data set, and W is the weighted adjacency matrix in the weighted adjacency matrix data set. Store the calculation result as a Laplacian matrix data set.

[0060] Step S324, solve the generalized eigenvalue problem. Use the Laplace matrix L in the Laplace matrix data set and the degree matrix D in the degree matrix data set to solve the equation L * φ_k = λ_k * D * φ_k. Store the obtained eigenvalues ​​and eigenvectors as an eigenvalue data set and an eigenvector data set, respectively.

[0061] Step S325, eigenvector selection. Sort the eigenvalues ​​in the eigenvalue data set, select the eigenvectors corresponding to the K smallest non-zero eigenvalues, and store the selected eigenvectors as the selected eigenvector data set.

[0062] Step S326: construct a low-dimensional embedding matrix. Combine the K eigenvectors in the selected eigenvector dataset into a matrix, with each column corresponding to one eigenvector. Store this matrix as an embedding matrix dataset.

[0063] Step S327: Calculate low-dimensional feature representation. Multiply each vertex in step S31 by the embedding matrix in the embedding matrix dataset to obtain the coordinates of each vertex in the low-dimensional space. Store the calculation result as a low-dimensional spatiotemporal feature dataset.

[0064] In this embodiment, first, the Laplace matrix of the graph is constructed. This method not only considers the connection relationship between data points, but also reflects the strength of spatiotemporal similarity through weights, effectively capturing complex spatiotemporal structures. Secondly, the generalized eigenvalue problem is solved. This method can extract the most important structural information in the data, and can better retain nonlinear spatiotemporal relationships than traditional principal component analysis. Again, a low-dimensional embedding matrix is ​​constructed by selecting the eigenvectors corresponding to the K smallest non-zero eigenvalues. This method achieves effective dimensionality reduction while retaining the key features of the data. Finally, each vertex is mapped to a low-dimensional space to obtain a new feature representation, which provides more compact and information-rich data for subsequent analysis. In the study of watershed vegetation, this atlas embedding method is of great significance. It can effectively handle the high dimensional and nonlinear characteristics in vegetation change data and extract the key features that best reflect vegetation dynamics. This method is particularly suitable for processing ecological data with complex spatiotemporal dependencies, such as long time series of vegetation indices or large-scale spatial distribution data. Through dimensionality reduction, this method not only reduces the redundancy and noise of the data, but also retains the main patterns and trends of vegetation changes. This is crucial for understanding the long-term response of vegetation to climate change and human activities. In addition, because atlas embedding takes into account the spatiotemporal relationship between data points, it can better capture the spatial continuity and temporal evolution characteristics of vegetation change, providing a new way to identify key ecological processes and driving factors.

[0065] Step S33: Dynamic system modeling Based on the low-dimensional spatiotemporal feature data set obtained in step S32, a nonlinear autoregressive exogenous model is constructed. First, the normalized vegetation index sequence is extracted from the low-dimensional spatiotemporal feature data set as the dependent variable, and other features are used as independent variables. Then, the lag orders ny and nu of the model are set. Next, an input matrix X is constructed, including the lagged dependent variable and the current and lagged independent variables. Then, a nonlinear function such as a kernel method or a neural network is used to fit the relationship between X and the current dependent variable y. The optimal model structure and parameters are determined by cross-validation. Finally, the trained model is used to fit and predict each spatiotemporal position to obtain a vegetation dynamic response model data set.

[0066] Specifically, step S33 is further as follows: Step S331: extract the normalized vegetation index sequence from the low-dimensional spatiotemporal feature data set obtained in step S327 as the dependent variable, and store it as the NDVI sequence data set. Take other features as independent variables and store them as the independent variable data set.

[0067] Step S332: Determine the optimal lag order ny (dependent variable) and nu (independent variable) by using an information criterion (such as AIC or BIC), and store the determined lag order as a lag order data set.

[0068] Step S333: Use the data in the NDVI sequence data set to construct a lagged dependent variable matrix, and use the data in the independent variable data set to construct the current and lagged independent variable matrices. Merge these matrices to obtain a complete input matrix. Store the constructed input matrix as the model input data set.

[0069] Step S334: compare the performance of multiple nonlinear functions (such as radial basis function, multilayer perceptron, etc.), select the most suitable nonlinear function type, and store the selected function type as a nonlinear function type data set.

[0070] Step S335: Use the model input data set and the NDVI sequence data set to train a nonlinear autoregressive exogenous model. Use cross validation to determine the optimal model parameters. Store the trained model parameters as a model parameter data set.

[0071] Step S336: Use the parameters in the model parameter data set to predict the test set data. Calculate the root mean square error, determination coefficient and other evaluation indicators between the prediction results and the actual NDVI values. Store the evaluation results as the model evaluation data set.

[0072] Step S337: Use the trained model to fit and predict each spatiotemporal position, and store the fitting and prediction results as a vegetation dynamic response model dataset.

[0073] In this embodiment, first, the optimal lag order is determined using information criteria (such as AIC or BIC). This method can scientifically select the time window of the model, taking into account the importance of historical information and avoiding overfitting. Secondly, an input matrix containing lagged dependent variables and current and lagged independent variables is constructed. This structure can simultaneously consider the influence of vegetation's own dynamics and external driving factors. Thirdly, the performance of various nonlinear functions (such as radial basis functions, multilayer perceptrons, etc.) is compared to select the most suitable nonlinear function type. This method can better capture the complex nonlinear characteristics of vegetation response. Finally, the optimal model parameters are determined by cross-validation, and a comprehensive model evaluation is performed, which ensures the generalization ability and prediction reliability of the model. In the study of watershed vegetation, this dynamic system modeling method has important application value. It can accurately describe the dynamic response process of vegetation to climate change and human activities, including short-term fluctuations, long-term trends, lag effects and cumulative effects. Compared with traditional static or linear models, the NARX model can better handle the nonlinear relationships and complex time dependencies that are prevalent in ecosystems. This is crucial for understanding and predicting the long-term evolution of vegetation under the background of climate change. Especially when studying the impact of extreme climate events on vegetation, this model can capture the nonlinear process and critical points of vegetation recovery. In addition, because the model takes into account the combined effects of multiple driving factors, it can provide a more comprehensive explanation of vegetation changes and help identify key ecological processes and environmental thresholds.

[0074] In summary, in step S3: First, a spatiotemporal weighted graph is constructed to integrate spatial and temporal information into a unified framework. This method can effectively capture the spatial dependence and temporal continuity of vegetation changes. Secondly, the graph embedding technology is applied to retain the key spatiotemporal features through dimensionality reduction, while reducing the redundancy and noise of the data, improving the efficiency and accuracy of subsequent modeling. Finally, a nonlinear autoregressive exogenous model (NARX) is used to construct a vegetation dynamic response model. This model can handle complex nonlinear relationships and long-term dependencies, and is more suitable for describing the dynamic response of vegetation to environmental changes than traditional linear models. This spatiotemporal dynamic modeling method is particularly suitable for the study of vegetation change at the basin scale because it can simultaneously consider spatial heterogeneity and temporal variability. In practical applications, this method can more accurately predict the response of vegetation to climate change and human activities, including hysteresis effects and cumulative effects. In addition, because the model takes into account spatiotemporal correlations, it can provide more reliable long-term predictions, which is crucial for formulating long-term ecological management strategies and assessing the potential impacts of climate change.

[0075] Step S4: separation and quantification of multi-factor influences The vegetation dynamic response model established in step S3 is used to separate and quantify the impacts of climate change and human activities on vegetation. First, the time-based causal entropy method is applied to identify the causal relationship between climate variables and human activities on the normalized vegetation index. Then, the improved Granger causal decomposition method is used to decompose the changes in the normalized vegetation index into the impact of climate change and human activities. Next, the generalized additive model is applied to quantify the nonlinear and threshold effects of each factor. Finally, the dynamic time warping algorithm is used to analyze the time lag effect of each influencing factor.

[0076] Step S41: Apply the time-based causal entropy method to the vegetation dynamic response model dataset obtained in step S33. First, calculate the transfer entropy between each pair of variables. Then, determine the significance threshold through permutation test. Next, construct a causal network based on the significant transfer entropy. Then, use conditional transfer entropy to eliminate indirect causal relationships. Finally, obtain the direct causal relationship network of climate variables and human activities on the normalized vegetation index, and store it as a causal relationship dataset.

[0077] Specifically, step S41 is further as follows: Step S411: extract time series data from the vegetation dynamic response model data set obtained in step S337, and store it as a time series data set.

[0078] Step S412: For each pair of variables in the time series data set, calculate the transfer entropy TE_X→Y from variable X to variable Y. Store the calculation result as a transfer entropy matrix data set.

[0079] Step S413: For each transfer entropy value in the transfer entropy matrix data set, generate a null distribution by randomly permuting the time series. Calculate the position of the actual transfer entropy value in the null distribution to obtain a p-value. Store the p-value as a significance matrix data set.

[0080] Step S414: Based on the p-value in the significance matrix data set, select significant transfer entropy relationships to construct an initial causal network, and store the network structure as an initial causal network data set.

[0081] Step S415: Calculate the conditional transfer entropy for each causal link in the initial causal network data set, taking into account the influence of other variables, and store the calculation results as a conditional transfer entropy matrix data set.

[0082] Step S416: compare the conditional transfer entropy in the conditional transfer entropy matrix data set with the original transfer entropy in the transfer entropy matrix data set to remove possible indirect causal relationships. Store the updated network structure as a refined causal network data set.

[0083] Step S417: Based on the refined causal network data set, a direct causal relationship network between climate variables and human activities and normalized difference vegetation index is constructed, and the final network is stored as a causal relationship data set.

[0084] In this embodiment, first, the transfer entropy is calculated. This method can quantify the transfer of information from one time series to another, which can better reflect the true causal relationship than the traditional correlation analysis. Secondly, a permutation test is performed to determine the significance threshold. This method evaluates the statistical significance of the observed causal relationship by generating a null distribution, effectively controlling the false discovery rate. Thirdly, the initial causal network is constructed and the conditional transfer entropy is calculated. This method can eliminate indirect causal relationships and identify direct causal links. Finally, the final causal network is generated, which provides a comprehensive view for understanding the causal structure in the system. In the study of watershed vegetation, this causal discovery method has important scientific significance. It can identify true causal relationships, not just correlations, from complex ecological data, which is crucial to understanding how climate change and human activities affect vegetation. Especially in the face of multiple potential driving factors, this method can distinguish between direct and indirect effects, helping researchers identify key ecological processes and environmental thresholds. In addition, since time dependence is taken into account, this method is particularly suitable for studying the lagged response and cumulative effects of vegetation to environmental changes.

[0085] Step S42: Based on the causal relationship dataset obtained in step S41, the improved Granger causal decomposition method is used to separate the impact. First, a complete model is constructed, including all climate variables and human activity indicators. Then, a restricted model containing only climate variables is constructed. Next, the prediction errors of the complete model and the restricted model are calculated. Then, based on the difference in prediction errors, the relative contributions of climate change and human activities are calculated. Finally, the decomposition results are applied to each spatiotemporal location to obtain a climate change impact dataset and a human activity impact dataset.

[0086] Specifically, step S42 is further as follows: Step S421: extract climate variables and human activity indicators from the causal relationship data set obtained in step S417, and store them as a climate variable data set and a human activity indicator data set, respectively.

[0087] Step S422: Use all variables in the climate variable dataset and the human activity index dataset to construct a complete regression model for predicting NDVI, and store the model parameters as a complete model parameter dataset.

[0088] Step S423: construct a restricted regression model for predicting NDVI using only the variables in the climate variable dataset, and store the model parameters as a restricted model parameter dataset.

[0089] Step S424: Use the parameters in the complete model parameter dataset and the restricted model parameter dataset to predict NDVI respectively, calculate the prediction errors of the two models, and store the results as a prediction error dataset.

[0090] Step S425: Calculate the relative contributions of climate change and human activities based on the complete model and restricted model prediction errors in the prediction error dataset, and store the calculation results as a relative contribution dataset.

[0091] Step S426: Apply the results in the relative contribution dataset to each spatiotemporal location to obtain the spatial distribution of climate change and human activity impacts. The decomposition results are stored as a climate change impact dataset and a human activity impact dataset, respectively.

[0092] In this embodiment, first, a complete model containing all climate variables and human activity indicators is constructed. This comprehensive model structure can take into account the combined effects of multiple factors. Secondly, by constructing a restricted model that only contains climate variables, an effective separation of the impacts of climate change and human activities is achieved. Thirdly, the relative contribution is calculated based on the difference in prediction errors. This method can not only quantify the impact intensity of each factor, but also evaluate their relative importance. Finally, the decomposition results are applied to each spatiotemporal location to obtain a high-resolution spatial distribution map, which provides important information for understanding the spatial heterogeneity of the impact. In the study of watershed vegetation, this impact separation method has great practical significance.

[0093] Step S43, applying the generalized additive model to the climate change impact data set and the human activity impact data set separated in step S42. First, a smoothing function, such as B-spline or thin plate spline, is set for each influencing factor. Then, the penalized likelihood method is used to estimate the parameters of each smoothing function. Next, the optimal smoothing parameter is determined by cross-validation. Then, the partial effect map of each factor is calculated. Finally, nonlinear relationships and threshold effects are identified based on the partial effect map to obtain a nonlinear effect data set.

[0094] Specifically, step S43 is further as follows: Step S431: extracting various influencing factor data from the climate change impact data set and the human activity impact data set obtained in step S426, and storing them as an influencing factor data set.

[0095] Step S432: For each factor in the influencing factor data set, try the B-spline and thin plate spline functions respectively, select the optimal smoothing function type through cross-validation, and store the selection result as a smoothing function type data set.

[0096] Step S433: Use the factors in the influencing factor data set as independent variables and NDVI as the dependent variable to establish an initial GAM model structure and store the model structure as a GAM structure data set.

[0097] Step S434: For the model in the GAM structure data set, use the penalty maximum likelihood method to estimate the parameters of each smoothing function, and store the estimated parameters as a smoothing function parameter data set.

[0098] Step S435: Use the generalized cross validation (GCV) or restricted maximum likelihood (REML) method to determine the optimal smoothing parameter for each smoothing function in the smoothing function parameter data set, and store the optimized smoothing parameter as the optimized smoothing parameter data set.

[0099] Step S436: Use the parameters in the optimized smoothing parameter data set to perform residual analysis, partial autocorrelation analysis, and other diagnostics on the model. The diagnostic results are stored as a model diagnostic data set.

[0100] Step S437: Calculate the partial effect diagram of each influencing factor based on the smoothing function parameter data set and the optimized smoothing parameter data set, and store the calculation result as the partial effect diagram data set.

[0101] Step S438: Analyze the curve shape in the partial effect diagram data set, identify significant nonlinear relationships and potential threshold points, and store the analysis results as a nonlinear effect data set.

[0102] In this embodiment, first, the optimal smoothing function is selected for each influencing factor. This method can better capture the nonlinear relationship between environmental factors and vegetation responses. Secondly, the penalized maximum likelihood method is used to estimate the smoothing function parameters. This method can strike a balance between model complexity and goodness of fit and avoid overfitting. Thirdly, the smoothing parameters are optimized by generalized cross validation (GCV) or restricted maximum likelihood (REML) methods to ensure the stability and generalization ability of the model. Finally, by calculating the partial effect diagram and analyzing the curve shape, significant nonlinear relationships and potential threshold points are identified, which provides key information for understanding the complex dynamics of the ecosystem. In the study of watershed vegetation, this nonlinear effect quantification method has important scientific value. It can reveal the complex patterns of the impact of environmental factors on vegetation, including non-monotonic relationships, saturation effects, and mutation points.

[0103] Step S44, use the dynamic time warping algorithm to analyze the time lag effect of each influencing factor obtained in steps S42 and S43. First, calculate the dynamic time warping distance matrix for each pair of time series. Then, find the optimal alignment path based on the distance matrix. Next, calculate the time difference on the alignment path to obtain the lag time distribution. Then, perform statistical analysis on the lag time distribution to determine the main lag mode. Finally, organize the analysis results into a time lag effect matrix and store it as a lag effect data set.

[0104] First, a time-based causal entropy method is used for causal discovery. This method can identify nonlinear causal relationships in complex systems and better reflect the true causal mechanism than traditional correlation analysis. Secondly, the improved Granger causal decomposition method is used for impact separation. This method not only considers the temporal relationship between variables, but also quantifies the relative contribution of different factors. Thirdly, the generalized additive model (GAM) is used to quantify nonlinear effects. This method can capture the complex nonlinear relationship and potential threshold effect of environmental factors on vegetation. Finally, the dynamic time warping (DTW) algorithm is used to analyze the time lag effect, which is crucial for understanding the delayed response of vegetation to environmental changes. In the study of watershed vegetation, this multi-factor impact separation and quantification method can solve the problem that has long plagued ecologists, that is, how to accurately distinguish the impact of climate change and human activities on vegetation. It can not only quantify the relative importance of each factor, but also identify the key impact threshold and time lag, which is of great guiding significance for the formulation of targeted ecological protection and restoration strategies. In addition, the flexibility of this method enables it to adapt to different types of watershed ecosystems, providing the possibility for cross-regional comparative studies.

[0105] Step S5: Based on the results of the previous steps, multi-scale impact assessment, scenario simulation, uncertainty quantification and comprehensive index construction are performed. First, the wavelet transform multi-resolution analysis method is used to evaluate the impact of climate change and human activities at different spatiotemporal scales. Then, a Monte Carlo simulation framework is constructed to generate multi-scenario predictions of future climate change and human activity impacts. Next, the Bayesian model averaging method is applied to quantify the uncertainty of the prediction. Finally, a comprehensive evaluation index system is constructed, including the impact intensity index, the response sensitivity index and the uncertainty index.

[0106] Step S51, applying the wavelet transform multi-resolution analysis method to the climate change impact dataset and the human activity impact dataset obtained in step S42. First, select a suitable wavelet basis function, such as Daubechies wavelet. Then, perform discrete wavelet transform on each time series. Next, reconstruct the signal at different decomposition levels to obtain impact sequences of different time scales. Then, perform statistical analysis on the reconstructed sequence, such as trend analysis and periodicity analysis. Finally, combine the analysis results of different time scales with spatial information to obtain a multi-scale impact assessment dataset.

[0107] Specifically, step S51 is further as follows: Step S511: extract time series data from the climate change impact data set and the human activity impact data set obtained in step S426, and store them as an impact time series data set.

[0108] Step S512: Compare the performance of different wavelet basis functions (such as Haar, Daubechies, Symlet, etc.) on the time series data set, select the most suitable wavelet basis function, and store the selection result as a wavelet basis function data set.

[0109] Step S513: Determine the maximum level of wavelet decomposition based on the length of the time series data set and the research objective, and store the determined number of levels as a decomposition level data set.

[0110] Step S514: Use the basis functions in the wavelet basis function data set and the number of levels in the decomposition level data set to perform discrete wavelet transform on each time series in the impact time series data set, and store the transform result as a wavelet coefficient data set.

[0111] Step S515: Use the coefficients in the wavelet coefficient data set to reconstruct the impact sequences of different time scales (such as inter-annual, inter-decadal, etc.), and store the reconstruction results as a multi-scale impact sequence data set.

[0112] Step S516: Perform trend analysis on each scale sequence in the multi-scale impact sequence data set, such as linear regression or Mann-Kendall test, and store the analysis results as a trend analysis data set.

[0113] Step S517: Perform periodic analysis on each scale sequence in the multi-scale impact sequence data set, such as Fourier analysis or wavelet power spectrum analysis, and store the analysis results as a periodic analysis data set.

[0114] Step S518: Combine the results of the trend analysis dataset and the periodicity analysis dataset to perform spatial cluster analysis on the study area, identify areas with similar multi-scale characteristics, and store the clustering results as a spatial clustering dataset.

[0115] Step S519: Integrate the results of the multi-scale impact sequence dataset, the trend analysis dataset, the periodicity analysis dataset, and the spatial clustering dataset to generate a detailed multi-scale impact assessment report, and store the report as a multi-scale impact assessment dataset.

[0116] In this embodiment, first, the most suitable wavelet basis function is selected, which ensures that the decomposition result can best capture the characteristics of the time series. Secondly, a discrete wavelet transform is performed to decompose the time series into components of different scales, which enables researchers to analyze short-term fluctuations and long-term trends at the same time. Thirdly, trend analysis and periodicity analysis are performed on sequences of different time scales, which helps to identify the main change patterns at different time scales. Finally, through spatial clustering analysis, the multi-scale analysis results are combined with spatial information to identify areas with similar multi-scale characteristics, which provides a new perspective for understanding the spatial pattern of vegetation change. In the study of watershed vegetation, this multi-scale impact assessment method has important theoretical and practical significance. It can simultaneously reveal the short-term fluctuations, medium-term cycles and long-term trends of vegetation changes, which is crucial for a comprehensive understanding of the response of ecosystems to environmental changes. For example, by distinguishing between interannual changes and long-term trends, researchers can more accurately assess the cumulative effects of climate change without being disturbed by short-term fluctuations. In addition, by identifying the dominant factors at different time scales, this method helps to explain complex ecological processes. For example, short-term vegetation changes may be mainly affected by weather conditions, while long-term trends may reflect more the impact of climate change or land use changes.

[0117] Step S52, using the vegetation dynamic response model dataset obtained in step S33 and the multi-scale impact assessment dataset of step S51, construct a Monte Carlo simulation framework. First, set possible scenarios of climate change and human activities based on historical data and future predictions. Then, add random disturbances to the input variables of each scenario to generate multiple sets of simulation inputs. Next, substitute the simulation inputs into the vegetation dynamic response model to obtain the corresponding vegetation change predictions. Then, repeat the simulation process multiple times to obtain the probability distribution of the prediction results. Finally, perform statistical analysis on the simulation results to obtain a multi-scenario prediction dataset.

[0118] Step S53, applying the Bayesian model averaging method to the multi-scenario prediction data set obtained in step S52. First, define a set of candidate models, including vegetation dynamic response models with different structures and parameters. Then, calculate the marginal likelihood of each model. Next, calculate the posterior probability of each model based on the marginal likelihood. Then, use the posterior probability as a weight to perform a weighted average on the prediction results of different models. Finally, calculate the variance of the weighted average prediction as a measure of uncertainty to obtain an uncertainty quantification data set.

[0119] Specifically, step S53 is further as follows: Step S531 : extracting model structure and parameter information from the vegetation dynamic response model data set obtained in step S337 , and storing it as a benchmark model data set.

[0120] Step S532: Based on the benchmark model data set, a series of candidate models are generated by changing the model structure (such as adding or deleting variables, changing the nonlinear function type, etc.) and the generated model set is stored as the candidate model data set.

[0121] Step S533: For each model in the candidate model data set, use the Laplace approximation or Monte Carlo integration method to calculate its marginal likelihood, and store the calculation results as a marginal likelihood data set.

[0122] Step S534: Calculate the posterior probability of each model based on the marginal likelihood values ​​in the marginal likelihood data set, and store the calculation results as a posterior probability data set.

[0123] Step S535: For each model in the candidate model dataset, prediction is performed using the multi-scenario prediction dataset obtained in step S52. The prediction results of all models are stored as a model prediction dataset.

[0124] Step S536: Use the posterior probabilities in the posterior probability data set as weights to perform weighted averaging on the prediction results in the model prediction data set, and store the average results as the BMA prediction data set.

[0125] Step S537: Based on the BMA prediction data set and the model prediction data set, calculate the variance of the prediction result as a measure of prediction uncertainty, and store the calculation result as a prediction variance data set.

[0126] Step S538: decompose the overall uncertainty in the prediction variance data set into model structure uncertainty, parameter uncertainty and scenario uncertainty, and store the decomposition result as an uncertainty decomposition data set.

[0127] Step S539: Integrate the results of the BMA prediction data set, the prediction variance data set, and the uncertainty decomposition data set to generate a detailed uncertainty quantification report, and store the report as an uncertainty quantification data set.

[0128] In this embodiment, first, a set of candidate models covering different structures and parameters is defined. This method takes into account the uncertainty of model selection. Secondly, the marginal likelihood and posterior probability of each model are calculated, which provides an objective method to evaluate the relative importance of different models. Thirdly, the posterior probability is used as a weight for model averaging. This method not only provides more reliable predictions, but also takes into account the differences between models. Finally, by decomposing the overall uncertainty into model structure uncertainty, parameter uncertainty and scenario uncertainty, it provides in-depth insights for a comprehensive understanding of the reliability of the prediction results. In the study of watershed vegetation, this uncertainty quantification method has great scientific and practical significance. It can provide more reliable and robust predictions of vegetation changes, while clearly indicating the uncertainty range of the prediction results. This is crucial for risk-based decision making, especially in the face of complex ecosystems and uncertain future climate scenarios.

[0129] Step S54, based on the multi-scale impact assessment data set of step S51, the multi-scenario prediction data set of step S52 and the uncertainty quantification data set of step S53, a comprehensive evaluation index system is constructed. First, the impact intensity index is calculated to quantify the impact of climate change and human activities on vegetation. Then, the response sensitivity index is calculated to indicate the sensitivity of vegetation to different influencing factors. Next, the uncertainty index is calculated to reflect the reliability of the prediction results. Then, these indices are standardized and weighted. Finally, the weighted indices are combined into an overall evaluation index to obtain a comprehensive evaluation index data set.

[0130] First, a multi-scale impact assessment was conducted using the wavelet transform multi-resolution analysis method, which can simultaneously analyze the change patterns at different time scales, thereby revealing the relationship between short-term fluctuations and long-term trends. Second, a Monte Carlo simulation framework was constructed for scenario simulation, which can generate multiple possible future scenarios and provide decision makers with a comprehensive risk assessment. Finally, the Bayesian model averaging (BMA) method was used for uncertainty quantification, which not only considers the uncertainty of the model structure, but also evaluates the uncertainty of parameter estimates and future scenarios. This comprehensive assessment and uncertainty analysis method has important practical significance in watershed vegetation management. It can help decision makers understand the dynamics of vegetation changes at different time scales, evaluate the effects of different management strategies under various possible scenarios, and quantify the reliability of prediction results. Especially in the face of the complex interactions between climate change and human activities, this method provides a scientific and comprehensive decision support tool. By constructing a comprehensive assessment indicator system, including impact intensity index, response sensitivity index and uncertainty index, this method provides a quantitative assessment standard for the sustainable management of watershed ecosystems, which helps to formulate more accurate and effective ecological protection policies.

[0131] In summary, in this embodiment: First, the acquisition and preprocessing of multi-source data provide a high-quality, spatially and temporally consistent data foundation for subsequent analysis. Secondly, the advanced feature extraction and fusion method takes into account multiple factors such as topography, land use changes, and changes in atmospheric composition, and comprehensively depicts the impact of environmental changes on vegetation. Spatiotemporal dynamic modeling innovatively integrates spatial and temporal information to accurately describe the dynamic response of vegetation to environmental changes. The multi-factor impact separation and quantification method not only distinguishes the impacts of climate change and human activities, but also quantifies the relative contribution and nonlinear effects of each factor. Finally, the comprehensive assessment and uncertainty analysis provides multi-scale, multi-scenario assessment results and quantifies the uncertainty of the prediction. This comprehensive and systematic approach significantly improves the ability to understand the dynamics of complex ecosystems.

[0132] Example 2: Obtain NDVI time series data (D_NDVI), vegetation functional type data (D_PFT) and vegetation coverage data (D_VCF) of the study area from the remote sensing database.

[0133] The precipitation (D_P), temperature (D_T) and relative humidity (D_RH) data of the study area were obtained from the meteorological database.

[0134] The digital elevation model (D_DEM), land use change data (D_LUC) and Concentration data( ).

[0135] All the data obtained above are processed for spatiotemporal consistency, including time series interpolation, spatial resolution resampling and projection transformation, to obtain a preprocessed data set (D_PRE) in a unified format.

[0136] Using the above temperature D_T and relative humidity D_RH, the improved Magnus formula is used to calculate the saturated water vapor pressure difference VPD (D_VPD). The formula is as follows: VPD = es*(1-RH / 100) Where es=0.611*exp((17.27*T) / (T+237.3)) T is temperature (°C) and RH is relative humidity (%).

[0137] Based on the digital elevation model D_DEM, the multi-scale terrain position index (MTPI) algorithm is used to extract the terrain feature D_TOPO. The MTPI calculation formula is: MTPI = Σ(wi * (Z0 - Zi)) / Σwi, where Z0 is the center pixel elevation, Zi is the neighborhood pixel elevation, and wi is the weight.

[0138] The improved Markov chain model is applied to the land use change data D_LUC to extract the land use change characteristics D_LUCC. The transition probability matrix P is calculated as follows: P = [pij], pij = nij / Σnij Where nij is the number of pixels that change from type i to type j.

[0139] Based on S14 Concentration data and NDVI time series data D_NDVI, using a nonlinear response function Quantification Fertilization effect : . Where a and b are fitting parameters, For reference concentration.

[0140] The preprocessed dataset D_PRE in a unified format is combined with the extracted features using the tensor decomposition method. After fusion, a high-dimensional feature tensor D_HFT is obtained.

[0141] Based on the high-dimensional feature tensor D_HFT, a spatiotemporal weighted graph G(V,E) is constructed, where V is the vertex set (representing spatial location) and E is the edge set (representing spatiotemporal relationship). Edge weight w ij The calculation is as follows: ij = exp(-d ij 2 / σ 2 ) * exp(-|t i -t j | / τ) where d ij is the spatial distance, t i and t j is the time point, σ and τ are adjustable parameters.

[0142] Apply the spectral clustering algorithm to the constructed spatiotemporal weighted graph G to obtain the low-dimensional spatiotemporal feature representation D_STFE. The eigenvector φ_k satisfies: L * φ_k = λ_k * D * φ_k where L is the Laplacian matrix of the graph, D is the degree matrix, and λ_k is the kth eigenvalue.

[0143] Based on the low-dimensional spatiotemporal feature representation D_STFE, a nonlinear autoregressive exogenous model (NARX) is used to construct a vegetation dynamic response model: y(t) = f(y(t-1), ..., y(t-ny), u(t-1), ..., u(t-nu)) + ε(t) where y is NDVI, u is exogenous input (climate factor), f is a nonlinear function, and ε is a noise term.

[0144] The time-based causal entropy (TE) method is applied to the dynamic model to identify the causal relationship between climate variables and human activities on NDVI. TE is calculated as follows: TE_Y→X = H(X|X-) - H(X|X-,Y-) where H is the conditional entropy, and X- and Y- represent the historical states of X and Y.

[0145] Based on the causal relationship, the improved Granger causal decomposition method is used to decompose the NDVI change into the impact of climate change (D_CCI) and the impact of human activities (D_HAI).

[0146] The generalized additive model GAM is applied to the above separated impacts (climate change impact (D_CCI) and human activity impact (D_HAI)) to quantify the nonlinear and threshold effects of each factor: NDVI = s1(X1) + s2(X2) + ... + sn(Xn)+ ε where si is the smoothing function, Xi is the influencing factor, and ε is the residual.

[0147] The dynamic time warping (DTW) algorithm is used to analyze the time lag effect of each influencing factor and obtain the lag effect moment (D_LEM).

[0148] Based on the analysis results, the wavelet transform multi-resolution analysis (WMRA) method was used to evaluate the impact of climate change and human activities at different temporal and spatial scales, and obtain multi-scale assessment results.

[0149] Using dynamic models and multi-scale assessment results, a Monte Carlo simulation framework is constructed to generate a multi-scenario prediction D_MSP of future climate change and the impact of human activities.

[0150] The Bayesian model averaging method is applied to the simulation results to quantify the uncertainty of the prediction. The BMA posterior distribution is calculated as follows: p(y|D) = Σ p(y|Mk,D) * p(Mk|D) where Mk is the model, D is the observed data, and p(Mk|D) is the model weight.

[0151] Based on the above results, a comprehensive evaluation index system was constructed, including the impact intensity index III, the response sensitivity index RSI and the uncertainty index UI.

[0152] Example 3: This example provides a method for quantifying the impact of climate change and human activities on grasslands in cross-border river basins in southwest China. Based on the normalized difference vegetation index NDVI, plant functional type PFT and vegetation continuous field VCF time series datasets, combined with improved vegetation coverage and coefficient of variation CV, a PCGP method for extracting grassland pixels with pure climate impact is constructed.

[0153] In order to extract grassland pixels that are as pure as possible, the present invention obtains pixels with a grassland coverage rate greater than 90% based on a commonly used grassland pixel extraction method.

[0154] In order to separate the pixels caused by climate, the sensitivity of CV growth change based on multi-year VCF dataset was used. Considering the important influence of VPD on vegetation, VPD was added as one of the independent variables to each pixel, and a multivariate regression model of NDVI with pure climate influence was established.

[0155] The trend analysis method was used to analyze the changing trends of NDVI and climate variables, and the multiple residual trend analysis method was applied to calculate the impact of human factors on grassland NDVI. Then, the partial correlation analysis method was used to eliminate the interference effects of other variables.

[0156] The evaluated model was applied to the effective grassland pixel VGP in the entire study area, and the NDVI (NDVI_c) inverted by the model was subtracted from the grassland NDVI observation data (NDVI_o) to obtain the NDVI value (NDVI_h) affected by human activities (Ha) to study the impact of climate change (Cc) and Ha on grassland vegetation in STB.

[0157] For each pixel, the present invention ignores areas with little or no grassland vegetation (annual average NDVI is less than 0.1); at the same time, pixels with NDVI>=0.1 and PFT>=75% are selected as VGP and used as valid NDVI_o pixels.

[0158] In order to clarify the NDVI pixels affected by climate and human activities, this paper proposes a new method to extract PCGP. Specifically: Based on VGP, pixels with PFT >= 90% are defined as pure grassland pixels (PGP); pixels with CV values ​​obtained from VCF data less than a certain threshold (CV < 19% in this study) are selected; pixels with CV < 19% overlapping with PGP are defined as PCGP, which is used to establish the NDVI model affected by climate. Among them, CV is the ratio of the standard deviation of the VCF time series of each pixel to the mean value, which is the most widely used indicator to measure the degree of fluctuation and is used to distinguish grasslands affected by climate and those not affected by climate.

[0159] Furthermore, according to the selected PGP, four climate variables are imported from the CMFD dataset: Precipitation, near-surface air temperature, surface air pressure, and near-surface air specific humidity; VPD is calculated from temperature, air pressure, and specific humidity; A multiple regression model between NDVI and climate variables (precipitation, temperature, and VPD) was established and evaluated using a 10-fold cross-validation method. Based on the model, the NDVI_c of each pixel was simulated. The difference between NDVI_c and NDVI_o is the residual (NDVI_h), which represents the response of the grassland to Ha, and the formula is as follows: , , Among them, NDVI_c, NDVI_h and NDVI_o refer to climate impact, human impact and observed grassland, respectively; P, Te and VPD refer to the average annual precipitation, temperature and VPD, respectively, in mm, °C and hPa; α, β and γ are the regression coefficients of multiple linear regression, and the average value of all simulations is taken; δ is the error term.

[0160] From the results, it can be seen that the effective grassland pixels in STB are mainly distributed in areas above 27°N, mainly concentrated in the middle and upper reaches of the three major international river basins. PGP mainly exists in areas with an altitude of about 4000-5200 meters. As we all know, the smaller the CV value, the more stable it is.

[0161] In PGP, the overall CV value is greater than 17%, and the CV value less than 19% is mainly distributed in the middle reaches of the Yarlung Zangbo River and the upper reaches of the Nujiang River and the Lancang River. The number of pixels is large and the distribution is balanced. The area with a CV value greater than 19% is mainly concentrated in the Yarlung Zangbo River, and the number of pixels decreases as the CV value increases.

[0162] Therefore, the present invention selects PGP with a CV value less than 19% for post-modeling.

[0163] The NDVIc obtained by the multivariate linear regression model based on PCGP (i.e., Cc-affected NDVI model) was closer to the observed value (NDVIo), with a CC of 0.79, a BIAS of −1.0%, and a RMSE of 0.10. Therefore, the model was reasonable and reliable and could be applied to VPG and distinguish the NDVI affected by Ha and Cc during the study period.

[0164] According to the research results, from 1982 to 1999, NDVI of grassland in 91.3% of the study area showed an upward trend, and the distribution of image pixels greater than 0.001 / a was more, especially in the Yarlung Zangbo River, followed by the Nujiang River. In the range of 0.001-0.002 / a, the trend range is very significant, and only scattered areas in the upper reaches of the Nujiang River have a slight downward trend. On the contrary, during the period of 2000-2015, 84.7% of the study area showed a downward trend, which may be closely related to the drought in 2009 / 2010 that destroyed vegetation growth in the southwest. The area with a clear increasing trend is located below the upper reaches of the Nujiang River (a small area at the end of the Yarlung Zangbo River).

[0165] In order to analyze the trend of NDVI in STB grassland, the present invention established a multiple linear regression model of precipitation, temperature, VPD and NDVI. During 1982-1999, the trend of NDVI change was greater than 0.001 / a, and more than 80% of the areas were significant.

[0166] In contrast, NDVI showed a significant downward trend over a large area from 2000 to 2015, especially in the middle of the Yarlung Zangbo River Basin. From 1982 to 1999, the NDVI changes in grasslands caused by Ha mainly showed a downward trend, ranging from 0.002 to 0.001 / a, with the largest area. Although NDVI showed a mainly downward trend from 2000 to 2015, the central area of ​​the Yarlung Zangbo River Basin showed a significant upward trend.

[0167] The study found that grassland showed an increasing trend from 1982 to 1999, with Cc and Ha (Ha) accounting for 73.15% and 5.74% of the increase in grassland NDVI, respectively. In addition, grassland degradation was only affected by Ha, which may be related to the construction of hydropower stations.

[0168] In contrast, grasslands in the Yarlung Zangbo River Basin showed a decreasing trend from 2000 to 2015, with Cc and Ha accounting for 33.30% and 25.00% of the grassland NDVI reduction area, respectively. Grasslands affected by Cc were most widely distributed in the Yarlung Zangbo River Basin, while grasslands affected by Ha were most prominent in the Lancang River Basin. In addition, grasslands driven by Cc and Ha (Cc-Ha) accounted for 21.30% of the total area, mainly concentrated in the Nujiang River Basin. Overall, in the MOST, Cc had a greater impact on grassland ecosystems than Ha.

[0169] In addition, the present invention also analyzes the relative contribution of Cc and Ha to grassland changes, as follows: During 1982-1999, Cc and Ha accounted for 80.87% and 19.13% of the total grassland change in STB, respectively. This indicates that grassland change in the region is mainly caused by Cc. In addition, 85.90% of grassland area in the Science and Technology Park was affected more by Cc than by Ha, and the trend was increasing. During 2000-2015, Cc and Ha accounted for 53.16% and 46.84% of the total grassland change in the Saudi Beach, respectively, indicating that grassland change in the region was mainly driven by Cc and Ha. In addition, the grassland area affected by Cc more than Ha accounted for 45.50% of the total grassland area in the region, and the trend was decreasing; this situation mainly occurred in the Yarlung Zangbo River Basin and the upper reaches of the Nujiang River Basin. At the same time, in STB, 40.20% of grassland area was affected more by Ha than by Cc, and the trend was decreasing; this situation mainly occurred in the Lancang River Basin.

[0170] The above results show that Cc is the main driving force of NDVI changes in STB grassland. In order to further explore the impact of climate variables on grassland NDVI changes, the present invention draws a spatial pattern map of driving forces based on the partial correlation coefficients of temperature, precipitation and VPD on NDVI and their identification criteria. For STB, the impact of climate variables on grassland is recovery from 1982 to 1999 and degradation from 2000 to 2015, which is consistent with the above results.

[0171] The area of ​​grassland driven by climate variables showed an increasing trend from 1982 to 1999. Grassland driven by temperature and precipitation (TP) was distributed throughout the Saudi Basin, accounting for 33.15%, followed by grassland driven by temperature (T), accounting for 29.44%.

[0172] In addition, grasslands driven by VPD are mainly distributed between 82°E and 86°E, concentrated in the Yarlung Zangbo River Basin; grasslands driven by precipitation (P) are mainly distributed around 89°E. In contrast, grasslands driven by climate variables showed a decreasing trend from 2000 to 2015. Grasslands driven by precipitation and VPD (PV) account for 20.56% of the STB area, followed by VPD, which accounts for 13.70%.

[0173] In general, in STB, grasslands were affected by multiple variables more than by a single climate variable: TP had the most significant effect on grassland recovery from 1982 to 1999, while PV had the most significant effect on grassland suppression from 2000 to 2015. In terms of the effects of single variables, climate variables restored grasslands from 1982 to 1999. The effects were followed by temperature, VPD, and precipitation. However, climate variables suppressed grasslands from 2000 to 2015, with the effects being followed by VPD, precipitation, and temperature.

[0174] From the above, it can be seen that the key step of the present invention is to extract the pixels caused by pure Cc. From the research results, it can be seen that the method proposed by the present invention is reasonable and effective, and can greatly improve the simulation results.

[0175] In addition, it should be noted that this paper only selected three climate variables (precipitation, temperature and VPD) to explore the climate sensitivity of grassland NDVI changes, which leaves room for uncertainty.

[0176] In addition, the present invention considers VPD as an important factor affecting grassland changes besides precipitation and temperature, and it can also be replaced by other influential climatic factors (such as nitrogen deposition or soil moisture); when evaluating the impact of climate on grassland normalized difference vegetation index, the number of climate components that should be included has not yet been determined, and sensitive variables and dominant variables have not yet been determined.

[0177] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solutions of the present invention may be modified or replaced by equivalents without departing from the spirit and scope of the technical solutions of the present invention, which should all be included in the scope of the claims of the present invention.

Claims

1. A method for quantifying the impact of climate change and human activities on watershed vegetation, characterized by: include, Obtain data from multiple data sources in the study area, perform spatiotemporal consistency processing, obtain preprocessed data sets in a unified format, and calculate the saturated vapor pressure difference; Extract features based on the preprocessed data set, and fuse the features with the preprocessed data set to construct a high-dimensional feature tensor; Based on the high-dimensional feature tensor, a spatiotemporal graph is constructed, and the dimension is reduced and a dynamic response model is established by using graph embedding technology to analyze the spatiotemporal dynamic characteristics of watershed vegetation; Causal network and regression decomposition methods were used to isolate and quantify the direct effects of climate and human activities on vegetation; Conduct a comprehensive assessment and perform uncertainty analysis.

2. A method for quantifying the impact of climate change and human activities on watershed vegetation as claimed in claim 1, characterized in that: The calculation of the saturated vapor pressure difference includes: Calculation of saturated water vapor pressure using the modified Magnus formula , the formula is: ; Wherein, T is the temperature data; Calculate the actual vapor pressure using the calculated saturated vapor pressure and relative humidity data , the formula is: ; Among them, RH is the relative humidity data; Calculating Saturation Vapor Pressure Deficit , the formula is: VPD = es - ea.

3. A method for quantifying the impact of climate change and human activities on watershed vegetation as claimed in claim 2, characterized in that: The extraction of features and fusion of the features with the processed data to construct a high-dimensional feature tensor includes the following steps: The MTPI algorithm is used to extract terrain features based on multi-scale windows; Applying the improved Markov chain model to extract the characteristics of land use change; Based on CO2 concentration data and vegetation index data, a nonlinear response function is used to quantify the CO2 fertilization effect. The processed data are fused with the extracted terrain features, land use change features, and CO2 fertilization effect features to obtain a high-dimensional feature tensor.

4. A method for quantifying the impact of climate change and human activities on watershed vegetation as claimed in claim 3, characterized in that: The method of extracting land use change characteristics using the improved Markov chain model comprises the following steps: Extracting land use change data from the preprocessed data set and storing it as a land use time series data set; For every two adjacent time points in the land use time series data set, the transition frequency between the land use types is calculated, the transition frequency is divided by the total frequency of the initial state to obtain a transition probability matrix P, and the calculated transition probability is stored as a transition probability matrix data set; For each spatial pixel, using the probability in the transition probability matrix data set, repeatedly simulate the land use change sequence within the study period to generate multiple change sequences, and store the generated change sequences as a land use change trajectory data set; Applying a sequence clustering algorithm to the land use change trajectory dataset, determining an optimal number of clusters, assigning a cluster label to each spatial pixel, and storing the clustering results as a land use change pattern dataset; Merging the land use change pattern dataset with the land use time series dataset, calculating the statistical characteristics of each change pattern, and storing the merged data as a land use change characteristic dataset; The land use time series dataset is a collection of extracted land use change data.

5. A method for quantifying the impact of climate change and human activities on watershed vegetation as claimed in claim 4, characterized in that: The method of quantifying the carbon dioxide fertilization effect using a nonlinear response function based on carbon dioxide concentration data and vegetation index data comprises the following steps: Extracting carbon dioxide concentration data and normalized vegetation index data from the preprocessed data set, and storing them as a carbon dioxide concentration data set and an NDVI data set, respectively; Construct a nonlinear response function model, where the nonlinear response function is: ， Among them, a and b are the parameters to be estimated. is the carbon dioxide concentration, is the reference value of carbon dioxide concentration; For each spatial pixel, a nonlinear least square method is used, with the logarithmic concentration ratio data set as an independent variable and the NDVI data set as a dependent variable, to fit the nonlinear response function, and the estimated a and b parameters are stored as a parameter estimation data set; Using the parameters in the parameter estimation dataset, calculate the predicted NDVI value for each spatiotemporal location; The parameters in the parameter estimation dataset and the logarithmic concentration ratio dataset are used to calculate the carbon dioxide fertilization effect at each spatiotemporal location, and the calculation results are stored as a carbon dioxide fertilization effect dataset.

6. A method for quantifying the impact of climate change and human activities on watershed vegetation as claimed in claim 5, characterized in that: The construction process of the high-dimensional feature tensor is as follows: Aligning the preprocessing dataset, the terrain feature dataset, the land use change feature dataset, and the carbon dioxide fertilization effect dataset in time and space dimensions, and storing them as an aligned feature dataset; reorganizing the aligned feature dataset into a three-dimensional tensor, and storing the constructed tensor as an original feature tensor dataset; Performing Tucker decomposition on the original feature tensor data set, initializing the core tensor and factor matrix, and iteratively optimizing the decomposition result using alternating least squares method until convergence or reaching a maximum number of iterations; Use the cross-validation method to determine the optimal decomposition rank and store it as a decomposition rank data set; Using the optimal rank in the decomposition rank dataset, performing a final Tucker decomposition on the original feature tensor dataset, and storing the obtained core tensor and factor matrix as a decomposition result dataset; The high-dimensional feature tensor is reconstructed using the core tensor and factor matrix in the decomposition result data set, and the reconstructed tensor is flattened into a two-dimensional matrix, in which each row represents a spatiotemporal position and each column represents a fused feature, and the two-dimensional matrix is ​​stored as a high-dimensional feature data set.

7. A method for quantifying the impact of climate change and human activities on watershed vegetation as claimed in claim 6, characterized in that: Constructing a spatiotemporal graph based on the high-dimensional feature tensor includes: Extract spatial coordinates and time information from the high-dimensional feature data set, and construct a spatiotemporal weighted graph G(V,E), where V is a vertex set and E is an edge set. Take each spatiotemporal position as a vertex of the spatiotemporal weighted graph and calculate the weight of the edge: ， in, is the spatial distance, and is the time point, σ and τ are adjustable parameters; The optimal σ and τ values ​​are determined by cross-validation, and a weighted adjacency matrix is ​​constructed based on the calculated weights to obtain a spatiotemporal weighted graph data structure.

8. A method for quantifying the impact of climate change and human activities on watershed vegetation as claimed in claim 7, characterized in that: The causal network and regression decomposition methods are used to separate and quantify the direct impacts of climate and human activities on vegetation, including: Use causal entropy method to construct causal network and quantify the direct causal relationship between climate and human activities; Separate and quantify the contribution of climate and human activities to vegetation through regression decomposition methods; The generalized additive model was used to quantify the nonlinear effects and thresholds of each influencing factor; The dynamic time warping algorithm is used to analyze the time lag effect and generate the lag effect matrix.

Images