Vegetation gross primary productivity monitoring method based on downscaling and multi-source remote sensing data

Abstract

The invention discloses a vegetation gross primary productivity monitoring method based on downscaling and multi-source remote sensing data, and belongs to the field of remote sensing image processing. In order to accurately monitor the GPP change of an engineering scale, the vegetation total primary productivity monitoring method comprises the following steps: acquiring first resolution multisource remote sensing data, downscaling by using CNN, inputting data containing DEM and LUCC in a layer, and dynamically adjusting background weight in a loss function through an LUCC binarization mask (a vegetation region is 1, and a non-vegetation region is 0); a CASA method is used for simulating GPP, a multi-year target area GPP data set is generated, and the maximum light energy utilization rate is determined by the vegetation type; analyzing GPP spatial and temporal change characteristics of the data set by using a Theil-Sen slope estimation method and the like, and extracting spatial distribution and a time change mode; extracting principal components of similar factors by using a principal component analysis method, calculating local correlation coefficients of the principal components and the GPP by using a partial correlation method, identifying dominant influence factors by using a contribution degree decomposition model, and extracting an environmental factor driven map. The method is applied to the field of ecological remote sensing monitoring.

Description

Technical Field

[0001] The present invention relates to the field of remote sensing information processing, and in particular to a method for monitoring gross primary productivity of vegetation based on downscaling and multi-source remote sensing data. Background Art

[0002] As hydropower development continues to advance, water conservancy and hydropower projects are delivering significant benefits in flood control and disaster reduction, clean energy supply, and ecological regulation. Hydropower has become the second largest conventional energy source. Statistics show that the 12 western provinces (autonomous regions, and municipalities) account for approximately 80% of the total hydropower reserves. This regional resource endowment has laid a solid foundation for cascade development of river basins and the construction of large-scale hydropower bases, effectively promoting the large-scale development and utilization of hydropower resources.

[0003] However, current hydropower development suffers from a lack of coordination between development planning and ecological and environmental protection. In some areas, overdevelopment has led to a series of negative effects, including ecosystem damage and ecological space compression, posing severe challenges to ecological protection. These issues have drawn significant attention from all sectors of society to the ecological impacts of hydropower development.

[0004] As a key component of global ecosystems, vegetation plays an indispensable role in climate regulation, shaping the surface environment, maintaining carbon balance, and regulating the water cycle. Gross primary productivity (GPP) refers to the total amount of organic carbon fixed by plants through photosynthesis at the ecosystem level. As the starting point of the terrestrial carbon cycle, GPP is a core indicator for measuring the dynamics of carbon sources and sinks and ecological regulation, and is of great significance for studying regional climate change and carbon balance. The dynamic evolution of GPP is driven not only by climate factors but also by human activities. The impacts of climate change and human activities on GPP exhibit dual characteristics. Therefore, continuously and accurately monitoring the spatiotemporal variations in GPP in terrestrial ecosystems under the influence of climate change and human activities will not only help us understand the dynamic evolution of ecosystems under the influence of climate change and human activities, and scientifically understand the intrinsic connections between ecosystems, human activities, and climate change, but also provide data support and scientific evidence for ecosystem management, laying a theoretical foundation for achieving sustainable development goals.

[0005] Monitoring plant carbon sequestration capacity can be traced back to the 18th century, when researchers focused on studying crop energy metabolism and ecological adaptation mechanisms. With the iterative upgrade of observation technology, the current technical paths for GPP monitoring have formed three categories: (1) ground-based measurement technology, (2) model deduction systems, and (3) canopy optical detection methods.

[0006] (1) For ground-based measurement technology, this technology has evolved from the early biomass harvesting method to a comprehensive observation system that can accurately extract vegetation carbon sequestration efficiency, obtain ecosystem dynamic parameters, and provide benchmark data for method validation. It mainly relies on biomass census and eddy covariance technology. Biomass census was the mainstream monitoring method in the 20th century. It is easy to operate but has problems such as long monitoring cycles. It is mostly used for local method parameter calibration. Eddy covariance technology is the core technology for carbon exchange research. Its advantage is non-destructive continuous monitoring. There are a large number of application cases, but there are bottlenecks in spatial scalability. It is mainly used for long-term ecological process analysis and as a benchmark for method validation.

[0007] (2) For model deduction systems, GPP monitoring in global or regional terrestrial ecosystems faces the challenge of lacking large-scale verification data. In recent years, a variety of methods have been applied. According to the simulation principle, the methods can be divided into four categories: (1) Statistical methods: Extraction methods based on the relationship between climate variables and biomass, such as the Miami method, have a simple structure but a single factor and high uncertainty. (2) Light energy utilization efficiency methods: Based on light energy utilization efficiency, there is a clear physiological basis and it is widely used, but the results may be inconsistent due to parameter differences. Common light energy utilization efficiency methods include VPRM, CE-LUE, MODE17, and the CASA method proposed in the existing technology 1. (3) Process-based methods: Comprehensive environmental factors, leaf level and ecosystem level, the mechanism is clear but requires a lot of data and complex parameterization. (4) Machine learning methods: Summarize relationships through data-driven, with high accuracy but black box characteristics, making it difficult to explain ecological mechanisms.

[0008] Prior art 1: Potter CS, Randerson JT, Field CB, et al., 1993. Terrestrial ecosystem production: A process model based on global satellite and surface data [J]. Global biogeochemical cycles, 7: 811-841.

[0009] Prior Art 2: Zhu Wenquan. 2005. Remote sensing estimation of net primary productivity of terrestrial ecosystem vegetation and its relationship with climate change[D]. Beijing Normal University.

[0010] As in prior art 2, the traditional CASA method sets the maximum light energy utilization rate to a fixed value of 0.389 g C·MJ⁻¹. The inventors found that this simplified treatment may underestimate the impact of vegetation type and habitat differences on light energy conversion efficiency.

[0011] (3) Regarding canopy optical detection methods, methods based on vegetation canopy reflectance use radiation characteristics to infer photosynthesis to monitor GPP, including sunlight-induced chlorophyll fluorescence (SIF) and near-infrared reflectance (NIRv). The former has a strong correlation with GPP, and the latter can generate a long-term global GPP dataset, but is limited in resolution and is only applicable to large scales. Each method has its own advantages and disadvantages: statistical methods are simple but have weak mechanisms, light energy utilization methods have a physiological basis but are greatly affected by the environment, process-based methods have clear mechanisms but high data requirements, machine learning methods have high accuracy but poor interpretability, and canopy reflectance methods have great potential but need to improve resolution. In research, it is necessary to combine regional vegetation characteristics, method parameters, and the balance between accuracy to monitor GPP.

[0012] In recent years, domestic and international research on the ecological and environmental impacts of water conservancy and hydropower projects has largely employed a method that combines satellite remote sensing data with measured vegetation data in the hydropower project areas, focusing on changes in land use, vegetation cover, and regional climate factors surrounding hydropower stations. However, existing vegetation remote sensing data has low spatial resolution and is discontinuous in time, making it difficult to accurately extract the specific ecological and environmental impacts of hydropower station construction at small scales, thus failing to accurately monitor the gross primary productivity of vegetation at these scales.

[0013] Furthermore, existing technologies have limitations in analyzing GPP impact characteristics. Traditional partial correlation analysis often relies on mathematical and statistical data processing methods (such as the three-factor control method). While these methods can eliminate the interference of specific variables, they struggle to extract spatial heterogeneity under the synergistic effects of multiple factors. Existing methods often prioritize climate as the core driving factor, while ignoring the impact of human activities on vegetation physiological processes. Furthermore, due to the limitations of graphical visualization technology, their analysis dimensions are often limited to three variables, resulting in systematic biases in the extraction of factor contributions.

[0014] Based on at least the above-mentioned status of the existing technology, at least one purpose of the present invention is to accurately monitor the impact of hydropower station construction on the temporal and spatial changes of GPP, and to make up for the various shortcomings of the relevant existing technology regarding the monitoring method of gross primary productivity of vegetation in small-scale engineering areas. Summary of the Invention

[0015] In order to alleviate or partially alleviate the above technical problems, the solution of the present invention is as follows:

[0016] A method for monitoring gross primary productivity of vegetation based on downscaling and multi-source remote sensing data includes the following steps:

[0017] Step S1, obtaining multi-source remote sensing data of a first resolution, and downscaling the multi-source remote sensing data using a CNN; wherein, in the process of downscaling the multi-source remote sensing data using the CNN, the input layer of the CNN includes at least DEM data and LUCC data, and the background area weight is dynamically adjusted in the loss function using a binary mask of the LUCC data, and in the binary mask of the LUCC data, the vegetation area is 1 and the non-vegetation area is 0;

[0018] Step S2, simulating vegetation GPP using the CASA method to generate a multi-year GPP dataset for the target area, wherein the maximum light energy utilization rate parameter in the CASA method is determined by the vegetation type;

[0019] Step S3: For the GPP dataset of the target area over many years, the temporal and spatial variation characteristics of GPP are quantitatively analyzed using the Theil-Sen slope estimation method, the Mann-Kendall non-parametric test method, and the coefficient of variation analysis method, and the spatial distribution and temporal variation pattern of GPP are extracted pixel by pixel;

[0020] In step S4, principal components of similar factors in environmental variables are extracted according to the principal component analysis method. Then, the local correlation coefficients of each principal component and the GPP are calculated pixel by pixel using the partial correlation method within the GPP range. Finally, the dominant influencing factors of each pixel are identified through the contribution decomposition model, and the environmental factor driving map at the engineering scale is extracted.

[0021] Furthermore, the method further includes step S0 of extracting the digital river network of the watershed, dividing the watershed area, and determining the watershed area controlled by the ridge line as the target area affected by the engineering construction.

[0022] Furthermore, the hydrological analysis module of ArcGIS was used to adopt the slope runoff simulation method and set different regional thresholds in the process of extracting the digital river network of the basin.

[0023] Furthermore, the extraction of the digital river network of the watershed and the division of the watershed area specifically include:

[0024] In ArcGIS, the Depression Fill tool of the Spatial Analyst module is called to determine the depression depth through iterative calculation and set the filling threshold. The depression elevation is raised to the lowest value of the adjacent unit to form a depression-free DEM.

[0025] The flow direction is determined by calculating the distance weight difference between the central grid and the neighboring cells, and the output is a raster data containing 1~255 flow direction codes;

[0026] Based on the no-sag DEM and flow direction data, the flow accumulation tool is used to count the number of upstream water catchment units of each grid to generate the flow accumulation raster;

[0027] A catchment area threshold is set to distinguish between river channels and non-river channels, and grids with accumulation greater than or equal to the threshold are identified as river networks.

[0028] Taking the valley as the pour point, the upstream watershed area is traced back in combination with the flow direction data, the Capture Pour Point tool is called to locate the outlet, and then the Watershed Tool is used to extract all the watershed grids in the basin to generate the watershed boundary extracted by vector.

[0029] Furthermore, the multi-source remote sensing data is unified to a first resolution by a bilinear interpolation method to obtain multi-source remote sensing data of the first resolution.

[0030] Furthermore, in the contribution decomposition model, the spatial distribution data of the interannual variation of GPP and environmental factors were first obtained through PCA and multi-factor partial correlation analysis;

[0031] Then, the contribution value estimation method is used to extract the contribution of the interannual variation of environmental factors to the interannual variation of GPP in the target area affected by engineering construction on a pixel-by-pixel basis.

[0032] Furthermore, the maximum patch index was selected from the patch level to represent the spatial proportion of the dominant landscape type, and the patch density was selected to reflect the degree of landscape fragmentation;

[0033] At the type level, the landscape shape index was selected to measure the complexity of patch morphology;

[0034] At the landscape level, the aggregation index was selected to evaluate the spatial connectivity of patches, and the Shannon diversity index and Shannon evenness index were used to extract the diversity characteristics of the ecosystem.

[0035] Furthermore, the multi-source remote sensing data includes terrain data and environmental factor data; the environmental variables are divided into four categories: temperature climate variables, moisture-related variables, vegetation characteristic variables and human activity representation variables.

[0036] Furthermore, the vegetation gross primary productivity monitoring method based on downscaling and multi-source remote sensing data is applied to an ecological remote sensing data processing system.

[0037] The technical method of the present invention has one or more of the following beneficial technical effects:

[0038] (1) This paper proposes a method based on digital river network extraction technology and the geomorphological characteristics of the study area. In an illustrative example, it systematically delineates vegetation areas affected by the construction of a cascade hydropower station on the Qinghai-Tibet Plateau. This method addresses the shortcomings of traditional hydropower station construction ecological research, which relies solely on artificially delineated large river basins or administrative divisions for regional selection, and effectively overcomes the uncertainty of results caused by the large target range.

[0039] (2) Although CNN performs well in downscaling tasks, existing downscaling methods often focus on a single vegetation index, temperature, or rainfall data, which may lead to the loss of small-scale features. This invention fuses LUCC data and DEM data, maintains the resolution of high-level feature maps, and integrates other types of data to extract a multi-factor coupled CNN downscaling method, breaking through the limitations of traditional single-variable downscaling methods and providing high-precision spatial data support for vegetation productivity monitoring in engineering construction areas. This invention verifies its applicability in downscaling long-term, multivariate remote sensing data and provides a reliable method for generating high-resolution environmental parameter products.

[0040] (3) Based on the inheritance of the classical method framework, the present invention focuses on the dynamic optimization of the maximum light energy utilization parameter, which can more accurately reflect the physiological characteristics of vegetation within the influence range of the hydropower station. It avoids the defect of the traditional CASA method that underestimates the impact of vegetation type and habitat differences on light energy conversion efficiency by simplifying the maximum light energy utilization parameter with a fixed value.

[0041] By downscaling remote sensing data with a temporally continuous 30-meter spatial resolution, combined with the CASA method, a more accurate engineering-scale GPP dataset can be extracted. This method systematically extracts the spatiotemporal dynamics of GPP and their impact characteristics before, during, and after the construction of a cascade of hydropower stations. This approach addresses the limitations of extracting the impacts of hydropower construction and other engineering-scale projects on small-scale vegetation and their impact characteristics.

[0042] Based on the digital river network extraction method and the downscaling processing method proposed in the present invention, and the dynamic optimization method for the maximum light energy utilization parameter in the CASA method, accurate data support and guarantee are provided for the accurate monitoring of GPP at the subsequent engineering scale.

[0043] (4) The four-dimensional environmental factor driving method extracted in this paper uses principal component analysis (PCA) to achieve dimensionality reduction, integration, and spatial expression of multi-source environmental factors, reclassifying environmental variables into four categories: temperature and climate variables, moisture-related variables, vegetation characteristic variables, and human activity-representing variables. On this basis, a three-level progressive analysis strategy is adopted to extract an engineering-scale environmental factor driving map. This method breaks through the spatial dimension limitations of traditional methods and can simultaneously analyze the spatiotemporal impacts of climate change and human activities on vegetation carbon sequestration.

[0044] In addition, other beneficial effects of the present invention will be mentioned in the specific embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1 This is a schematic diagram of the impact range of a certain cascade hydropower station;

[0046] Figure 2This is a diagram of the steps for obtaining the catchment area of ​​a cascade hydropower station basin;

[0047] Figure 3 1 is a schematic flow chart of the CASA method of the present invention;

[0048] Figure 4 1 is a diagram showing the accuracy comparison between the improved CASA method of the present invention and other existing technologies;

[0049] Figure 5 It is the result of CNN downscaling of multi-source remote sensing data in a certain embodiment of the present invention;

[0050] Figure 6 This is the spatial distribution characteristic map of total primary productivity;

[0051] Figure 7 This is the spatial distribution map of GPP changes based on the difference analysis of the mean GPP values ​​in the three years before and after construction;

[0052] Figure 8 is the spatial distribution of the interannual change rate of total primary productivity during the study period;

[0053] Figure 9 It is the spatial distribution map of the significance test of total primary productivity during the study;

[0054] Figure 10 is the spatial distribution map of the interannual variability of total primary productivity during the study period;

[0055] Figure 11 This is the result of the analysis of the factors affecting the spatial distribution of GPP;

[0056] Figure 12 This is a flow chart of the method for monitoring gross primary productivity of vegetation based on downscaling and multi-source remote sensing data of the present invention. DETAILED DESCRIPTION

[0057] To make the objectives, technical methods, and advantages of the present invention more clear, the technical methods of the present invention will be clearly and completely described below in conjunction with the accompanying drawings. Obviously, the embodiments described are part of the embodiments of the present invention, not all of them. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts should be within the scope of protection of the present invention.

[0058] In the description of the present invention, unless otherwise specified, “ / ” indicates that the objects associated with each other are in an “or” relationship. For example, A / B can represent A or B. “And / or” in the present invention is only a description of the association relationship of associated objects, indicating that there can be three relationships. For example, A and / or B can represent: A exists alone, A and B exist at the same time, and B exists alone. A and B can be singular or plural.

[0059] In the embodiments of the present invention, words such as "exemplarily" and "for example" are used to indicate examples, illustrations, or explanations. Any embodiment or design method described as "exemplarily" or "for example" in the embodiments of the present invention should not be construed as being preferred or advantageous over other embodiments or design methods. Rather, the use of words such as "exemplarily" and "for example" is intended to present the relevant concepts in a concrete manner to facilitate understanding.

[0060] The term "determining" in the present invention encompasses a variety of actions, and "determining" may include calculating, computing, processing, deriving, investigating, searching (e.g., via searching in a table, a database, or another data structure), ascertaining, etc. Also, "determining" may include receiving (such as receiving information), accessing (such as accessing data in a memory), etc. Furthermore, "determining" may include resolving, selecting, choosing, establishing, and other similar actions.

[0061] List of terms:

[0062] Multi-source remote sensing data: This includes, but is not limited to, terrain data (e.g., DEM data, LUCC data), environmental factor data (e.g., LST data, ET data), and other remote sensing data from various sources. Preferably, this data may also include future data (e.g., CMIP6 data).

[0063] Gross Primary Productivity (GPP) refers to the total amount of organic carbon fixed by plants through photosynthesis at the ecosystem level.

[0064] Land Use and Cover Change (LUCC): A key variable in land-atmosphere interactions, it strongly impacts Earth's ecological environment by altering surface properties. This research involves assessing land cover change; modeling and forecasting global land use and cover change; the linkages between the drivers of land use and cover change at global, regional, and local scales; and data development activities and information systems. For example, the CN Multi-Period Land Use Remote Sensing Monitoring Dataset (CNLUCC) uses a two-level classification system: a first-level classification system of six categories, primarily based on land resources and their use attributes, including cultivated land, forest land, grassland, water area, construction land, and unused land; and a second-level classification system of 23 types, primarily based on the natural attributes of land resources.

[0065] First resolution: refers to a resolution of remote sensing data. For example, it may be 30m×30m resolution. It may also be other suitable high resolutions. The present invention is not limited thereto.

[0066] Convolutional Neural Network (CNN): It is a neural network method that automatically learns low-, medium-, and high-level abstract features from raw pixels through multi-level feature extraction, effectively reducing the need for manual feature engineering.

[0067] As an example, the present invention takes the vegetation GPP in a certain cascade hydropower station construction area on the Qinghai-Tibet Plateau as the research object, and extracts the main processing steps covered by a specific embodiment of the method of the present invention. Those skilled in the art know that the steps involved in the method are also applicable to the processing of various types of remote sensing data about other research areas or objects. The present invention is not limited to this and is not limited to this.

[0068] The method for monitoring gross primary productivity of vegetation based on downscaling and multi-source remote sensing data, described in this invention, can be applied to ecological remote sensing data processing systems. These remote sensing image processing systems, or remote sensing data processing systems, are primarily used to preprocess, analyze, interpret, and apply remote sensing image data acquired from platforms such as satellites and aircraft. Their core functions include data correction, enhancement, classification, target detection, and information extraction. They are typically targeted at specific industries (such as meteorology, agriculture, geology, and ecology) and integrate industry-specific workflows.

[0069] Regarding the sources of multi-source remote sensing data for the research object, this invention utilizes three types of data: terrain data, environmental factor data, and future data. Terrain data and environmental factor data are primarily used to define the impact area of ​​a cascade hydropower station and study the spatiotemporal trends and impacts of GPP. Future data, on the other hand, is used to predict future spatiotemporal trends and their impact characteristics. As an example, Table 1 shows the three different data sources and their main characteristics.

[0070] Table 1: Sources of three different types of data in the present invention

[0071]

[0072] The Digital Elevation Method (DEM) data is derived from the global 30-meter resolution DEM data released by the NASA Shuttle Radar Topography Mission (SRTM). This dataset, acquired using interferometric synthetic aperture radar (InSAR) technology, covers the global land area from 60°N to 56°S. Its spatial resolution is 30 meters, with a vertical accuracy of ±16 meters absolute and ±10 meters relative. The data dates back to 2000 and is accessible through Google Earth Engine (GEE). Preprocessing, including stitching, cropping, reprojection, and quality control, was performed on this platform for terrain analysis, river mesh extraction, and CNN downscaling.

[0073] In an illustrative example, the time span from 2005 to 2022 can be selected, accessed through the GEE platform, and pre-processed by stitching, cropping, reprojection, and quality control to determine the light energy utilization parameters in the CASA method and CNN downscaling study.

[0074] The land use data are derived from the CLCD dataset (Table 2), which was released by Wuhan University in 2021.

[0075] Table 2: Classification of land use types in the CLCD dataset

[0076]

[0077] The monthly mean temperature data (LST) comes from the MYD11A1.061 dataset, which is released by the NASA EOSDIS Land Processes Distributed Active Archive Center (LP DAAC) and provides daily land surface temperature (LST) and emissivity information per pixel on a global scale.

[0078] The actual evapotranspiration (ET) and potential evapotranspiration (PET) evapotranspiration data are derived from the MOD16A2.061 dataset, which is a global land evapotranspiration data product obtained by the MODIS sensor onboard the Terra satellite of NASA's Earth Observing System (EOS). It belongs to the 6.1 version of the dataset (ArchiveSet 61).

[0079] The Downward Shortwave Radiation (DSR) data comes from MCD18A1.062, a dataset provided by NASA that combines observations from the Terra and Aqua satellites. It is designed to meet the needs of terrestrial methods and applications for high temporal and spatial resolution radiation data.

[0080] In addition, the fraction of photosynthetically absorbed radiation (FPAR) data is derived from the MOD15A2H.061 dataset, which is released by the NASA Earth Observing System (EOS) Land Processes Distributed Active Archive Center (LP DAAC).

[0081] The vegetation index data used in this paper include the Normalized Difference Vegetation Index (NDVI), the Enhanced Vegetation Index (EVI), and the Leaf Area Index (LAI). These indices reflect the ability of plants to absorb and store atmospheric carbon. The NDVI and EVI data are derived from the MODIS Land Standard Product MOD13Q1.061 (spatial resolution: 1 km, temporal resolution: 16 days), provided by NASA. The LAI data are derived from the MODIS Main Product MOD15A2H.061 (spatial resolution: 1 km, temporal resolution: 8 days).

[0082] In one embodiment of the present invention, the vegetation GPP data may be OD17A2H.061 GPP data and 500 m resolution GPP data of the Qinghai-Tibet Plateau.

[0083] In addition, the regional CMIP6 multi-model fusion gross primary productivity (GPP) dataset, released in 2022 by the Nanjing University of Information Science and Technology, covers the period 1850 to 2100, encompassing both the historical period (1850-2014) and future scenarios (2015-2100). The future data are extracted based on four shared socioeconomic pathways (SSP126, SSP245, SSP370, and SSP585) with a spatial resolution of 0.25°. This dataset utilizes deep learning methods, using the GLASS GPP data as a benchmark, and integrates the outputs of 23 CMIP6 models.

[0084] Climate data can be derived from the core scenario model of the Sixth Coupled Model Intercomparison Project (CMIP6), covering four representative shared socioeconomic pathways (SSPs) scenarios: SSP126, SSP245, SSP370, and SSP585 (Table 3).

[0085] Table 3: Four SSP scenarios and their differences

[0086]

[0087] For example, the future climate variables used can be the CMIP6 daily-scale pre, tmax, and tmin data (Table 4) output from the ACCESS-CM2 core scenario pathway for the period 2023–2100. These data were preprocessed in R Studio, including format conversion, cropping, reprojection, and quality control. The daily rainfall for each scenario was summed to obtain the monthly rainfall. The daily maximum temperature was calculated using the maximum method to obtain the monthly maximum temperature, and the daily minimum temperature was calculated using the minimum method to obtain the monthly minimum temperature. The monthly climate variable data for each scenario were then averaged and output.

[0088] Table 4: Description of future climate data in the CMIP6 dataset

[0089]

[0090] On the other hand, for the above data set, the multi-source remote sensing data processing method of the present invention includes the following steps:

[0091] Optionally, in step S0, the digital river network of the watershed is extracted, the watershed area is divided, and the watershed area controlled by the ridge line is determined as the target area affected by the engineering construction.

[0092] Furthermore, in the process of extracting the digital river network of the basin, the catchment area controlled by the ridge line is regarded as the target area affected by the construction of a hydropower station.

[0093] Defining the impact of hydropower station construction on regional vegetation and its scope (target area) is a key step in ecological impact assessment. Previous studies have shown that hydropower station operation leads to secondary succession of vegetation communities by altering the hydrological regime (fluctuations in upstream and downstream water levels and changes in evaporation patterns) and the soil matrix (saturation of soil in upstream submerged areas and increased erosion in downstream exposed areas). NDVI time-series analysis shows that within the 10 km buffer zone surrounding the Lancang River cascade development, vegetation cover exhibits spatiotemporal heterogeneity, with significant decreases in NDVI values ​​during the construction period and gradual increases during the recovery period. 63% of the area shows an ecological improvement trend, which captures the phased nature of the impact of hydropower projects on vegetation. The natural barrier effect of the river valley topography (average slope >25°) significantly limits the horizontal spread of vegetation impacts. This is supported by the spatial scale of the 80-meter rock mass deformation impact range, as demonstrated by monitoring in the Jinping Hydropower Station project area. In high-altitude ecologically fragile areas such as the Qinghai-Tibet Plateau, the coupling effects of low temperature, strong radiation environment and engineering disturbance will amplify the ecological response. For example, the reverse succession of vegetation that occurred during the construction of power stations in the middle reaches of the Yarlung Zangbo River requires the use of cold-resistant pioneer species for ecological restoration.

[0094] With respect to the exemplary embodiments of the present invention, taking into account the topographic constraints and ecological sensitivity of the unique alpine river valley areas of the Qinghai-Tibet Plateau, the inventors believe that the definition of the impact range of hydropower projects must integrate topography, climate, and ecology as a composite factor.

[0095] Based on this, the method for dividing the impact range of hydropower stations based on hydrogeological units proposed in this invention realizes the multi-dimensional coupling of terrain factors (slope, elevation), hydrological elements (runoff path, soil moisture) and ecological responses by extracting river networks and water catchment areas in the basin.

[0096] Specifically, the hydrological analysis module of ArcGIS can be used, and the slope runoff simulation method can be adopted to set different regional thresholds in the river network extraction process. Different thresholds lead to differences in the density of the water network. When the threshold is large, the cumulative number of grids required for water flow is large, and the river network is more difficult to form; when the threshold is small, the number of grids required for water flow is small, and the river network is easier to form. When setting the convergence accumulation threshold, the actual situation of the Yarlung Zangbo River system should be fully considered, and the threshold should be reasonably selected through observation and comparison with topographic maps, thereby ensuring the accuracy and practical applicability of the runoff simulation. When the upstream catchment area threshold is determined to be 5000 grid units, it can be seen from the overlay analysis with the first-level river network that the river system of a cascade hydropower station basin extracted by the present invention is basically consistent with the actual river network, which can meet the needs of hydrological simulation and water conservancy and hydropower development.

[0097] Topographic analysis results show that the terrain of a cascade hydropower station basin runs northwest to southeast, with lower elevations in the center and higher elevations on both sides, resulting in a steeper river channel slope. With the exception of peaks in the northwest, which exceed 5 km in elevation, the valley and southeast peaks range from 3 to 4 km in elevation. Digital river network extraction and analysis revealed that the spatial distribution of ridgelines and valley topography significantly influence the formation of watersheds. Specifically, the watersheds on either side of the valley exhibit significant spatial differentiation.

[0098] Figure 1 This is a schematic diagram of the impact area of ​​a cascade hydropower station. Within the study area, the distance from the ridgeline to the central axis of the river valley typically ranges from 2 km to 4 km. Due to the significant influence of the ridge watershed, the edges of the catchment area exhibit an irregular distribution that aligns closely with the ridgeline, forming independent hydrological units. This, in turn, results in a distinct "symmetrical segmentation" of the catchment area's spatial distribution. Because water is a key factor controlling vegetation growth, the high mountainous terrain on both sides weakens the hydropower station's ability to regulate climatic factors.

[0099] In one embodiment, the present invention considers the catchment area controlled by a ridgeline as the impact zone for the construction of a cascade of hydropower stations and serves as the target area. Compared to traditional buffer zone delineation methods, this catchment-based target area division not only considers direct impacts from hydropower station construction (such as inundation areas and construction sites), but also encompasses the secondary impacts of local microclimate changes on vegetation caused by microtopographic changes.

[0100] Furthermore, in order to extract the digital river network within the hydropower station construction basin, the DEM data can be preprocessed to determine the flow direction and the confluence accumulation calculation.

[0101] Figure 2 This is a diagram showing the steps for obtaining the catchment area of ​​a cascade hydropower station. To obtain the catchment area, a hydrological analysis method based on the digital elevation model (DEM) is used, combined with spatial analysis tools on the ArcGIS platform. The specific steps are as follows:

[0102] 1) Filling depressions

[0103] Raw DEM data often contain data depressions (due to instrument errors or truly closed depressions), which require filling to eliminate spurious terrain interference. Using the Fill tool in the Spatial Analyst module in ArcGIS, the depth of the depression is determined through iterative calculations. A fill threshold is set, raising the depression elevation to the lowest value of the adjacent cells to create a depression-free DEM. This step ensures continuous water flow and avoids interruptions in runoff simulation.

[0104] 2) Flow direction determination

[0105] The D8 single-direction flow algorithm simulates flow paths: water in each grid cell flows only in the direction of the steepest slope (i.e., the largest elevation difference) among its eight adjacent cells. This method determines flow direction by calculating the weighted distance difference (the ratio of elevation difference to grid spacing) between the central grid cell and its neighboring cells. The output is a raster containing flow direction codes ranging from 1 to 255. While this algorithm ignores multidirectional flow diversion and hydrological parameters (such as rainfall), it is computationally efficient and effectively reflects the dominant influence of topography on runoff.

[0106] 3) Calculation of cumulative flow

[0107] Based on the non-depression DEM and flow direction data, the Flow Accumulation tool was used to count the number of upstream catchment cells in each raster, generating a runoff accumulation raster. Accumulation values ​​indicate the intensity of runoff convergence: high values ​​correspond to areas with developed river channels, while zero values ​​represent watersheds. This step provides a basis for river network extraction.

[0108] 4) River network extraction

[0109] A catchment area threshold (5000, for example) is set to distinguish between river channels and non-river channels: Grids with a cumulative volume ≥ the threshold are considered river networks. The threshold selection requires a balance between river network density and topographic detail, typically determined through critical value analysis. After extracting the target raster using the Raster Calculator, the Raster River Vector Extraction tool is used to generate vector river networks, and the geometry is optimized through smoothing.

[0110] 5) Watershed generation

[0111] Using the valley of a cascade hydropower station as the pour point, we combined flow direction data to trace back to the upstream catchment area. We used the Capture Pour Point tool to precisely locate the outlet. Then, we used the Watershed tool to extract all the catchment rasters in the hydropower station's watershed, generating a vector-extracted catchment boundary. This method can clearly define the spatial extent of the watershed and provide key hydrological boundaries for water conservancy project planning.

[0112] Step S1, obtaining multi-source remote sensing data of a first resolution, and downscaling the multi-source remote sensing data through CNN; wherein, in the process of downscaling the multi-source remote sensing data using CNN, the input layer of CNN includes at least DEM data and LUCC data, and the background area weight is dynamically adjusted in the loss function through the binary mask of the LUCC data, wherein the vegetation area in the binary mask of the LUCC data is 1 and the non-vegetation area is 0.

[0113] Specifically, step S1 may include the following specific steps:

[0114] Multi-source remote sensing data were unified to a 30 m resolution using bilinear interpolation, and a binary LUCC mask (LUCC_Value) was established (with vegetation areas set to 1 and non-vegetation areas set to 0). This effectively eliminated the interference of background values ​​on weight calculation. Finally, the input layer of the CNN integrated four data channels, including DEM data, LUCC data, and their derived features.

[0115] For example, CNN can adopt the following three-layer convolution structure:

[0116] (a) Input layer (Conv2d): Channel number 4→64, kernel size 3×3, followed by adaptive average pooling (avgpool2d) for dimensionality reduction and preservation of critical space;

[0117] (b) Feature hidden layer (Conv2d): Channel number 64→32, kernel size 3×3, Dropout2d (dropout rate 0.2) is used to suppress overfitting;

[0118] (c) Output layer (Conv2d): Channel number 32→1, kernel size 1×1, resolution restoration achieved through upsampling, and 500 iterations of training using the adaptive minimum error method, with each iteration slicing the annual data.

[0119] To address the problem of edge pixel distortion, this paper introduces the LUCC binary mask LUCC_Value, which can dynamically adjust the background area weight in the loss function. The training process uses the backpropagation algorithm to update the convolution kernel parameters and minimize the mean square error through gradient descent. With the coefficient of determination (R²) and root mean square error (RMSE) as the core indicators, the multi-year mean comparison method is used to verify the spatiotemporal consistency of the downscaling results.

[0120]

[0121]

[0122] Where, represents the predicted value; represents the average value; represents the measured value; n represents the number of samples; R 2 It is a key parameter to measure the efficiency of the method, and its value range is between 0 and 1. A higher R 2 The value indicates that the CNN method has better performance. 2 When it is equal to 1, it indicates perfect fitting. RMSE is the root mean square error of the CNN method. A smaller RMSE indicates a higher method effectiveness.

[0123] Although the existing technology 3 has extracted the CNN-based terrain factor downscaling method, the root mean square error (RMSE) in some local areas is significantly lower than that of traditional methods (such as COSMO-CLM), verifying the effectiveness of CNN in geoscience data processing.

[0124] Existing technology 3: Quan Kai. 2021. Research on terrain factor downscaling method based on convolutional neural network[D]. Northwest Agriculture and Forestry University.

[0125] While CNNs excel at downscaling, they still have potential limitations. For example, the downsampling operation of traditional CNNs can lead to the loss of small-scale features. However, this method effectively alleviates this problem by maintaining the resolution of high-level feature maps and combining multiple sources of data (such as land use types).

[0126] Alternatively, dynamic modeling of multi-temporal information (such as the temporal fusion strategy of PCNN-GRU) can be used to enhance the ability to focus on features in key areas, thereby further improving the downscaling accuracy under complex surface conditions.

[0127] Alternatively, the attention mechanism can be introduced to enhance the ability to focus on features in key areas, thereby further improving the downscaling accuracy under complex surface conditions.

[0128] The CNN method of the present invention can break through the limitations of traditional univariate downscaling methods by fusing terrain factors (slope, elevation, etc.) with land use data, and provide high-precision spatial data support for vegetation productivity analysis in engineering construction areas.

[0129] In addition, the design of the CNN network architecture of the present invention takes into account both computational efficiency and feature retention requirements. The pooling-unpooling mechanism achieves a balance between feature compression and reconstruction, while the Dropout strategy enhances the adaptability of the method to irregular image structures.

[0130] The proposed CNN downscaling method, based on multi-source remote sensing data, demonstrated excellent performance. Linear regression analysis of all downscaled data with the original resolution data showed R² values ​​exceeding 0.9, and a highly significant correlation (P < 0.001) between the simulated data and the original data, demonstrating the strong explanatory power of the proposed method. Specifically, the RMSEs of the linear regressions of the downscaled and undownscaled data were: PAR 4.4 MJ·m⁻², FPAR 1.03%, ET 2.6 mm, PET 7.32 mm, and LST 0.93°C. For vegetation parameters, the RMSEs for NDVI, EVI, and LAI were 0.03, 0.01, and 0.45, respectively.

[0131] Based on this, the present invention uses multi-source remote sensing data after CNN downscaling for CASA method simulation and influencing factor analysis, which will have higher reliability.

[0132] Step S2: using the CASA method to simulate vegetation GPP, generating a multi-year GPP dataset for the target area, wherein the maximum light energy utilization parameter in the CASA method is determined by the vegetation type.

[0133] For example, the multi-year GPP dataset may be a dataset from 2005 to 2022, and further, may be data of the first resolution.

[0134] In the present invention, the maximum light energy utilization rate parameter in the CASA method is no longer a fixed value commonly used in the traditional CASA method.

[0135] Figure 3 The figure is a flow chart of the CASA method of the present invention. The CASA method proposed in prior art 1 has two outstanding advantages: (1) a simple parameter system, which only requires remote sensing data, meteorological data and vegetation index to drive it; (2) strong spatial scalability, which can achieve seamless simulation from site to large scale. After nearly three decades of development, this method has been widely used in fields such as global carbon cycle assessment and regional ecological monitoring. The entire document is incorporated into this application by reference. However, the traditional CASA method sets the maximum light energy utilization rate to a fixed value of 0.389 g C·MJ⁻¹. This simplified treatment may underestimate the impact of vegetation type and habitat differences on light energy conversion efficiency.

[0136] GPP is influenced by both the amount of photosynthetically active radiation absorbed by vegetation and the efficiency of light energy utilization. The amount of photosynthetically active radiation absorbed by vegetation is closely related to its absorption ratio and the total solar radiation, while the efficiency of light energy utilization is related to factors such as the maximum light energy utilization rate, the temperature stress coefficient, and the water stress coefficient.

[0137] Based on the assumption that vegetation has different maximum light energy utilization rates, the present invention combines vegetation classification data and the maximum light energy utilization rates of vegetation of corresponding classifications to dynamically optimize the maximum light energy utilization rate parameters. The CASA method can be expressed as follows:

[0138]

[0139] Where GPP represents the total carbon fixed by vegetation per unit area through photosynthesis (g C·m⁻²), The APAR is the total photosynthetically active radiation absorbed by the vegetation canopy (MJ m²), and is:

[0140] APAR=FPAR×SOL×0.5;

[0141] In the formula, SOL represents the total solar radiation (unit: MJ·m⁻²), FPAR is the ratio of photosynthetically active radiation absorbed by vegetation, and the constant 0.5 represents the ratio of photosynthetically active radiation to total solar radiation. In addition,

[0142]

[0143] Where, is the maximum light energy utilization rate of vegetation (g C·MJ⁻¹). The maximum light energy utilization rate of vegetation here is not a fixed value, but is determined according to the vegetation classification data. and is the temperature stress influence coefficient (temperature stress influence coefficient 1, temperature stress influence coefficient 2 respectively), is the water stress impact coefficient; and

[0144]

[0145] Where, Reflects the ability of vegetation to limit photosynthesis under high and low temperature environments. The optimum temperature for vegetation photosynthesis (degrees Celsius), usually the monthly average temperature when the NDVI value is the highest. When the monthly average temperature is less than or equal to -10℃, The value is 0.

[0146]

[0147] Where, The light energy utilization rate decreases when the ambient temperature changes from the optimum temperature to high or low temperature. is the monthly average temperature (℃). When the monthly average temperature is 13℃ lower or 10℃ higher than the optimum temperature, The value is half of the optimum temperature.

[0148]

[0149] Where, Reflects the effect of vegetation water conditions on light energy utilization, E is the actual evapotranspiration of the region (mm), is the regional potential evapotranspiration (mm).

[0150] For example, the performance evaluation indicators of the CASA method include the coefficient of determination (R²) and the root mean square error (RMSE). The correlation between the simulated GPP data and the existing GPP data is tested through linear regression analysis. The main reference dataset that can be used is the MODIS GPP dataset.

[0151] For example, the monthly mean temperature data (LST) is used to calculate the temperature stress index in the CASA model after preprocessing such as splicing, cropping, reprojection and quality control.

[0152] For example, actual evapotranspiration (ET) and potential evapotranspiration (PET) evapotranspiration data are used to calculate the water stress index in the CASA model after preprocessing such as stitching, cropping, reprojection and quality control.

[0153] For example, by using downlink shortwave radiation (DSR) data and preprocessing such as stitching, cropping, reprojection, and quality control, the photosynthetic active radiation (PAR) parameter required by the CASA model can be calculated by multiplying the DSR by 0.5.

[0154] For example, land use data is used, and after pre-processing such as stitching, cropping, reprojection and quality control, to determine the light energy utilization parameters and CNN downscaling process in the CASA model.

[0155] For example, the photosynthetic absorption ratio (FPAR) data is used and pre-processed by stitching, cropping, reprojection and quality control, and then used as the FPAR parameter in the CASA model.

[0156] For example, MODIS NDVI, EVI, and LAI data were used, and after multiple preprocessing steps such as stitching, cropping, reprojection, and quality control, these data were used as auxiliary parameters and influencing factors in the CASA method.

[0157] For example, the prediction accuracy of the CASA method can be verified using MODIS GPP data and 500 m resolution GPP data of the Qinghai-Tibet Plateau.

[0158] The above data can be obtained from at least the GEE platform, the National Tibetan Plateau Data Center, etc.

[0159] This paper improves the CASA method and integrates multi-source remote sensing environmental factors (including surface temperature, moisture, and landform type) to extract a high-precision method that takes into account human interference. Method validation shows that the linear regression coefficient R² of CASA-GPP, MODIS-GPP, and He-GPP data is greater than 0.85. The reason for this is that the present invention uses the dynamic maximum light energy utilization rate (ε g ) parameterization method, replacing the static vegetation type lookup table method used in MODIS-GPP, to extract the impact of hydropower development on land use processes. The results validate the effectiveness of GPP spatial accuracy in areas with significant human activities, such as hydropower interference.

[0160] Figure 4 The figure is a comparison of the accuracy between the improved CASA method of the present invention and other existing technologies. Figure 4 Part (a) is the comparison with MODIS-GPP. Figure 4 Part (b) is a comparison with He-GPP (prior art 4).

[0161] Prior art 4: He S, Zhang Y, Ma N, et al., 2022. A daily and 500 mcoupled evapotranspiration and gross primary production product across CNduring 2000–2020 [J]. Earth System Science Data Discussions, 2022: 1-42.

[0162] Results show that the improved CASA method, when used to extract GPP within the influence area of ​​a cascade hydropower station (at a spatial resolution of 30 m), exhibits good correlation with both MODIS-GPP and He-GPP linear regression validation results. Specifically, the CASA GPP obtained by this method has an R² of 0.913 and an RMSE of 90.09 g C m⁻² yr⁻¹ relative to MODIS GPP. Accuracy results demonstrate that the improved CASA method can accurately simulate the spatiotemporal variations in GPP within the hydropower station's influence area.

[0163] In step S3, for the GPP dataset of the target area over many years, the temporal and spatial variation characteristics of GPP are quantitatively analyzed using the Theil-Sen slope estimation method, the Mann-Kendall non-parametric test method, and the coefficient of variation analysis method, and the spatial distribution and temporal variation pattern of GPP are extracted pixel by pixel.

[0164] The Theil-Sen slope estimation method, a robust nonparametric statistical method, effectively extracts the direction and rate of interannual variation in GPP by calculating the median slope of each pixel's time series data. Its advantages are that it has no specific requirements for data distribution and is highly robust to measurement errors and outliers. The core of this method is to capture trend characteristics by calculating the median slope value of all pixel pairs within the study period.

[0165]

[0166] In mathematical expressions, represents the result of median operation on the slope sequence of n(n-1) / 2 inter-annual pixel combinations, and its numerical characteristics reflect the rate of change of the grid unit during the study period. Specifically, when When the value is positive, it indicates that the spatial unit presents an increasing trend; when When it is a negative value, it indicates a decreasing trend. and Respectively represent the first and The grid cell value corresponding to each observation year.

[0167] The complementary MK nonparametric test method determines the statistical significance of the trend at a significance level of 0.05 by standardizing the statistic Z value. This method does not rely on the normal distribution assumption and can overcome the interference of missing values ​​and outliers in time series.

[0168] The synergistic application of the Theil-Sen slope estimation method and the MK nonparametric test rule forms an analytical chain of "slope calculation-significance test": the Theil-Sen method provides trend strength extraction indicators, while the MK test gives it the rigor of statistical inference. This combined strategy maintains the adaptability of nonparametric methods and enhances the interpretability of trend analysis.

[0169] The Coefficient of Variation (CV) analysis method is used to extract the time series fluctuation characteristics of vegetation GPP. This indicator effectively characterizes the interannual variability of vegetation productivity through the relative ratio of the standard deviation to the mean, eliminating the influence of dimensional differences on the analysis results. When the CV value approaches zero, it indicates that vegetation productivity is highly stable. As the CV value increases, the volatility of the GPP series increases significantly.

[0170] The Theil-Sen slope estimation method, MK nonparametric test method, and CV analysis method are all common knowledge in the art, and more specific implementation details thereof will not be repeated in the present invention.

[0171] In step S4, principal components of similar factors are extracted using the principal component analysis method, and then the local correlation coefficient between each principal component and the GPP is calculated pixel by pixel using the partial correlation method. Finally, the dominant influencing factors of each pixel are identified through the contribution decomposition model, and the environmental factor driving map at the engineering scale is extracted.

[0172] The purpose of step S4 is to explore the dominant influencing factors of interannual variation of vegetation GPP.

[0173] Current research on the dynamics of GPP under the influence of hydropower projects has largely focused on low-altitude regions, while systematic studies of sensitive high-altitude ecological zones, particularly the Qinghai-Tibet Plateau, remain relatively scarce. Traditional methods for analyzing the association between environmental factors and GPP often employ bivariate correlation analysis using the Pearson correlation coefficient and the Spearman rank correlation coefficient. The Pearson correlation coefficient is suitable for continuous variables with linear relationships, while the Spearman rank correlation coefficient is less restrictive on data distribution and can capture monotonic trends.

[0174] However, the environmental changes caused by hydropower projects are characterized by multi-factor coupling. Conventional correlation coefficients are difficult to eliminate collinearity interference, which can easily lead to false correlations or underestimate the true strength of the association. Partial correlation analysis, however, controls for the statistical interference of other variables and accurately analyzes the "net correlation" between a single influencing factor and GPP.

[0175] This method has been validated in ecosystem carbon flux research (Prior Art 5) and effectively eliminates the impact of cross-talk from multiple environmental factors on analytical results. By extracting a multi-scale environmental factor dataset and combining it with a data processing framework for GPP data, this method facilitates the analysis of the pathways by which hydropower station construction impacts the carbon cycle within the unique geographical context of the Qinghai-Tibet Plateau.

[0176] Prior Art 5: Sun Hong, Fang Guofei, Ruan Linlin, et al. 2022. Spatiotemporal patterns and driving factors of carbon and water fluxes in semi-arid regions of Asia[J]. Acta Ecologica Sinica, 42(12):4742-4757.

[0177] To address the limitations of existing research in analyzing GPP impact characteristics, this paper proposes and applies an improved multi-factor spatial partial correlation analysis method. Traditional partial correlation analysis often relies on mathematical statistics (such as the three-factor control method). While this method can eliminate the interference of specific variables, it is difficult to extract spatial heterogeneity characteristics under the synergistic effects of multiple factors.

[0178] Existing methods often take climate factors as the core driving factors, while ignoring the impact of human activities on vegetation physiological processes. Moreover, due to the limitations of graphic visualization technology, their analysis dimensions are often limited to three variables, resulting in systematic deviations in the quantification of the contribution of influencing factors.

[0179] The four-dimensional driving force analysis framework extracted by this paper uses principal component analysis (PCA) to achieve dimensionality reduction, integration, and spatial expression of multi-source environmental factors. Specifically, this paper reclassifies environmental variables into four categories: temperature climate variables (LST corresponds to Temp), water-related variables (ET and PET correspond to Pre), vegetation characteristics (NDVI, EVI, and LAI correspond to VEG), and human activity indicators (Landscape Pattern Index corresponds to LSI).

[0180] On this basis, the present invention adopts a three-level progressive parsing strategy:

[0181] (a) First, principal components of similar factors are extracted through PCA to eliminate multicollinearity between variables;

[0182] (b) Secondly, the partial correlation method is used to calculate the local correlation coefficient between each principal component and GPP pixel by pixel;

[0183] (c) Finally, a contribution decomposition model is used to identify the dominant influencing factors for each pixel and extract a map of environmental factors driving the project-scale. This method transcends the spatial dimensionality limitations of traditional analysis and can simultaneously analyze the spatiotemporal impacts of climate change and human activities on vegetation carbon sequestration.

[0184] Furthermore, the present invention uses the Landscape Pattern Index (LPI) based on LUCC data to quantitatively analyze the evolution of spatial heterogeneity in the hydropower station construction impact zone. As a spatial metric characterizing landscape structure, the Landscape Pattern Index can capture the effects of human activities on landscape pattern across three dimensions: patch morphology, spatial configuration, and ecosystem function.

[0185] In the index selection process, in order to avoid indicator redundancy and highlight ecological significance, a hierarchical screening strategy was adopted: first, the maximum patch index (LPI) was selected from the patch level to characterize the spatial proportion of dominant landscape types, and the patch density (PD) reflected the degree of landscape fragmentation; secondly, the landscape shape index (LSI) was selected at the type level to measure the complexity of patch morphology; finally, the aggregation index (AI) was selected at the landscape level to evaluate the spatial connectivity of patches, and the Shannon diversity index (SHDI) and Shannon evenness index (SHEI) were used to jointly extract the diversity characteristics of the ecosystem.

[0186] Specifically, LPI identifies the dominant landscape type in a region by the proportion of the largest patch area (0-100%), and the PD value range is positively correlated with the degree of fragmentation; LSI characterizes the degree to which the morphology deviates from regular geometry by the ratio of patch perimeter to area; AI (0-100%) reflects the aggregation characteristics of patches of the same type, and high-value areas indicate stable ecological units; SHDI and SHEI analyze landscape heterogeneity from the perspectives of richness and balance, respectively, among which SHEI (0-1) has a significant negative correlation with dominance.

[0187] The Largest Patch Index (LPI) is used to characterize the spatial distribution and relative dominance of dominant landscape types within a study area. The index is measured on a continuous, closed scale from 0 to 100%, and its numerical characteristics clearly correlate with landscape pattern: a positive increase in the LPI reflects an increase in the spatial proportion of large-scale patches within the landscape matrix, while a decrease in the LPI indicates a decrease in the spatial dominance of large-scale patches within the study area.

[0188]

[0189] Where LPI represents the maximum plaque index, a max It represents the area of ​​the largest patch in a certain landscape type, and A represents the total area of ​​the landscape.

[0190] Patch density (PD) is an important parameter for quantifying landscape spatial pattern. It is defined as the spatial distribution density of specific landscape types within a study area. This index effectively characterizes the spatial heterogeneity and fragmentation of landscape systems: an increase in PD indicates increased spatial dispersion of landscape units within the study area, reflecting a higher level of landscape fragmentation; conversely, a decrease in PD indicates a trend toward spatial aggregation of landscape units and enhanced landscape connectivity.

[0191]

[0192] Where PD represents plaque density, N represents the total number of plaques, and A represents the total plaque area.

[0193] The Landscape Shape Index (LSI) is an important quantitative indicator of landscape spatial morphology, primarily used to assess the degree of deviation between the geometry of landscape patches and a standard square. This index effectively reflects the boundary complexity and spatial configuration characteristics of specific landscape types: a high LSI value indicates that patch edges are becoming more complex and the spatial structure is significantly irregular; conversely, a low LSI value indicates that the patch geometry approximates a regular square and the spatial structure is becoming more simplified.

[0194]

[0195] Where LSI represents the landscape shape index, E represents the total length of the patch boundary, and A represents the total area of ​​the patch.

[0196] The Aggregation Index (AI) is a key metric for quantifying landscape spatial configuration characteristics, used to assess the degree of spatial clustering and spatial associations among landscape patches. The index, measured on a closed scale from 0 to 100, exhibits a significant spatial correlation with landscape pattern: a decrease in the AI ​​value reflects an increase in the spatial clustering of landscape patches, accompanied by a decrease in landscape connectivity. Conversely, an increase in the AI ​​value indicates a trend toward a more homogeneous distribution of landscape units, with a significant increase in the degree of spatial aggregation among patches.

[0197]

[0198] Where AI represents the aggregation index, and gii represents the number of links between similar patches in the i-th type of landscape.

[0199] The Shannon's Diversity Index (SHDI) effectively characterizes the spatial heterogeneity of landscape systems. Changes in the index's values ​​are significantly correlated with land-use patterns: an increase in the SHDI indicates a simultaneous increase in the diversity of land-use types and the degree of spatial fragmentation within the study area. Numerous empirical studies have shown a significant positive correlation between landscape diversity and biodiversity, with this correlation exhibiting a typical normal distribution. Therefore, the SHDI not only reflects landscape pattern characteristics but also provides an important quantitative basis for assessing regional biodiversity levels.

[0200]

[0201] Where SHDI represents the Shannon diversity index, m represents the total number of patch types, and the Pi value represents the ratio of the area of ​​patch type i to the total landscape area.

[0202] The Shannon's Evenness Index (SHEI) is an important indicator for assessing the spatial distribution characteristics of landscapes, used to quantify the diversity and uniformity of landscape types. This index, which spans the closed interval [0,1], together with the Shannon Diversity Index, forms a comprehensive indicator system for landscape diversity analysis. Notably, the SHEI exhibits a significant negative correlation with the landscape dominance index: when the SHEI approaches 1, it indicates that a specific patch type dominates the landscape, exhibits significant spatial imbalance, and exhibits relatively low landscape diversity. Conversely, when the SHEI approaches 0, it indicates that patch types are evenly distributed across space, and landscape diversity remains high.

[0203]

[0204] Where SHEI represents the Shannon Evenness Index, m represents the total number of patch types, and the Pi value represents the ratio of the area of ​​patch type i to the total landscape area.

[0205] Principal component analysis is a statistical dimensionality reduction method based on the covariance matrix. Its core lies in converting the original variables into linearly independent principal components through orthogonal transformation, thereby reducing dimensional complexity while retaining the maximum data variation. Specifically, the PCA method extracts an orthogonal basis to linearly combine the original variables, maximizing the variance in the direction of the first principal component. Subsequent principal components then capture the remaining variation under orthogonal constraints, ultimately achieving data simplification by screening principal components with significant variance contributions. In this process, the principal component loading (i.e., the linear correlation coefficient between the original variable and the principal component) reflects the relative contribution of each variable to the formation of the principal component, as shown in the following formula:

[0206]

[0207] Where, express dimensional vector; express dimensional principal components; Indicates the principal components and The linear correlation coefficient of the variables, that is, the loading.

[0208] In order to further analyze the interaction of multiple variables, the present invention introduces multi-order partial correlation coefficient analysis. This method quantifies the net correlation between target variables by controlling the influence of a specific set of variables. When the number of control variables is three, the calculated partial correlation coefficient can effectively isolate the interference of the other three types of environmental factors and separately evaluate the effect of the remaining variables on GPP. The range of the correlation coefficient is [-1, 1]. The closer its absolute value is to 1, the more significant the linear correlation between the variables, and the sign represents the positive or negative direction of the effect. By systematically comparing the partial correlation coefficients under different control conditions, the dominant environmental factors that drive the interannual variation of NPP can be identified. The calculation formula of the multi-order partial correlation coefficient is as follows:

[0209]

[0210] In the formula, the partial correlation coefficient Represents the control variable 、 and When the linear effect of is eliminated, the research variable and The value range of this statistic is defined in the closed interval [-1,1], and its numerical characteristics have clear statistical significance: Indicates that there is a positive correlation between the two variables. Indicates that there is a negative correlation between the two variables, and the absolute value of the correlation coefficient is positively correlated with the strength of the association between the variables. The closer the value is to 0, the stronger the independence of the variables. The closer the value is to 1 or -1, it indicates that there is a significant linear dependence between the variables.

[0211] The GPP contribution estimation method first uses PCA and multi-factor partial correlation analysis to obtain the spatial distribution data of the interannual variation of GPP and environmental factors. Next, the contribution value estimation method is used to extract the contribution of the interannual variation of environmental factors to the interannual variation of GPP within the influence area of ​​the hydropower station on a pixel-by-pixel basis. The formula for the contribution estimation is as follows:

[0212]

[0213]

[0214] Where, For the The contribution of the interannual variation of GPP (or environmental factors) of each pixel to the interannual variation of GPP (or environmental factors); For the Pixels in Interannual variation of GPP (or environmental factors) at each moment; is the year of the data (for example, 2005 ≤ t ≤ 2022); It represents the interannual variation of GPP within the influence area of ​​the hydropower station.

[0215] Finally, the present invention uses a specific embodiment of some intermediate results or data involved in the processing of multi-source remote sensing data to demonstrate more details or effects of the present invention.

[0216] Figure 5 This is the result of CNN downscaling of multi-source remote sensing data in one embodiment of the present invention. In this embodiment, based on a convolutional neural network downscaling model, the present invention generated a monthly 30-meter spatial resolution dataset of atmospheric variables (LST, ET, PET, PAR, FPAR) and vegetation indices (NDVI, EVI, LAI) from 2005 to 2022. Compared to the original data, the downscaled data, while largely maintaining the original range, significantly enhances the representation of spatial detail, resulting in finer pixel structure and clearer topographic features and gradients.

[0217] Atmospheric variables exhibit significant spatial heterogeneity: PAR (Parity of the Earth's atmosphere) shows a gradient increasing from northwest to southeast, with monthly mean values ​​ranging from 173.41 to 206.11 MJ·m². FPAR (Four-Parity of the Earth's atmosphere) remains high in the upstream hillsides (monthly mean 7.4%-32.13%), while significantly decreasing in the downstream valleys. ET (Emergency of the Earth's atmosphere) is low in the upstream (monthly mean 15.26-39.49 mm) but generally higher in the downstream. PET (Potential Temperature of the Earth's atmosphere) is low in the northwest's high-elevation regions and higher in the lower elevations (monthly mean 50.37-116.22 mm). LST (Long-term Temperature Stress) is generally higher downstream than upstream, with an annual mean temperature difference of 15.62°C (5.49-21.11°C). In terms of vegetation indices, NDVI, EVI, and LAI all show a vertical differentiation characteristic of "low in river valleys - high on mountain slopes - low on mountain tops" (the annual average values ​​of NDVI range from 0.11 to 0.55, EVI is 0.05 to 0.28, and LAI is 1.33 to 8.32).

[0218] Figure 6 This is a spatial distribution map of gross primary productivity, specifically from 2005 to 2022. GPP within the hydropower station's impact area exhibits significant spatial heterogeneity. During the construction period from 2005 to 2022, mean GPP values ​​in the region ranged from 105.35 to 552.34 g C m⁻² yr⁻¹, with an overall average of 377.01 ± 79.12 g C m⁻² yr⁻¹, and an average annual total GPP of 0.0758 Tg C. The spatial distribution shows that GPP values ​​are closely related to altitude. Specifically, GPP values ​​in river valleys are generally lower than those on moderately elevated slopes, while the highest peaks in the northwest exhibit the lowest GPP values. Low GPP values ​​occur at the hydropower station's cutoff point, in the upstream valley, and in the downstream valley with low water levels. High GPP values ​​are mostly found on upstream slopes, some downstream slopes, and on low-elevation peaks.

[0219] Figure 7 This is a spatial distribution map of GPP changes, based on a difference analysis of the mean GPP values ​​before and after construction over the three-year period. This difference analysis reveals that GPP changes ranged from -165.51 to 163.28 g C m⁻² yr⁻¹, with an average change of 27.67 ± 34.47 g C m⁻² yr⁻¹. Significant decreases were primarily concentrated in the river valleys surrounding the hydropower station and upstream. Increases were observed on some hillsides and in downstream areas, while changes were relatively stable in other areas.

[0220] Figure 8This is the spatial distribution of the interannual rate of change of gross primary productivity during the study period. Based on the construction cycle characteristics of cascade hydropower stations, this example divides the period from 2005 to 2022 into different stages of project implementation for analysis. Considering the overlap between these stages during actual construction, this example divides the 18-year observation period into three phases based on the specific project progress: pre-construction (2005-2008), construction (2008-2019), and post-construction (2019-2022). Hydropower Station A was completed in 2014, and Hydropower Station B was under construction from 2015 to 2019. Due to the lag in remote sensing data updates, the post-construction phase only covers three years of data, from 2019 to 2022, making it difficult to independently analyze its spatiotemporal trends. To this end, this example will comprehensively analyze the spatiotemporal trends of GPP under the influence of cascade hydropower station construction in the following four time periods: (1) the entire study period (2005-2022) reflects the overall evolution; (2) the construction period (2008-2019) captures the direct impact of the project; (3) the impact of the project construction on the original ecology is extracted by combining the pre-construction and construction period (2005-2019); and (4) the construction period and post-construction period (2008-2022) observe the continued effect of the project.

[0221] Figure 9 This is the spatial distribution of gross primary productivity (GPP) significance tests during the study. Changes in GPP within the hydropower station's area of ​​influence during the entire study period, from 2005 to 2022, exhibited distinct spatiotemporal variations. Overall, 89.45% of the area showed a positive growth trend, while 10.55% experienced negative growth. Specifically, 41.2% of the area showed significant positive growth, while 2.52% showed significant negative growth. The interannual rate of change ranged from -3.79 to 6.56 g C m⁻² yr⁻¹, with an overall average growth rate of 2.45 g C m⁻² yr⁻¹. Negative growth areas were primarily concentrated in the upstream valley of the hydropower station's cutoff point, while positive growth areas were distributed along the upstream slopes and most of the downstream area.

[0222] Figure 10 This is a spatial distribution of the interannual variability in gross primary productivity during the study period. Based on the spatial distribution pattern of interannual variability in GPP, vegetation GPP fluctuations were low in most areas affected by hydropower construction between 2005 and 2022. Areas of low and relatively low fluctuation accounted for 0.3% and 56.92% of the total area, respectively (a total of 57.22%), while areas of moderate, relatively high, and very high fluctuation accounted for 38.11%, 3.71%, and 0.96% of the total area (a total of 42.78%), with CV values ​​ranging from 0 to 0.54. Areas of high fluctuation were primarily located near hydropower cutoffs, in the upstream river valley, and in some downstream areas, and their distribution was similar to that of areas with significant positive and negative growth.

[0223] Figure 11 This figure, based on the results of an analysis of factors influencing the spatial distribution of GPP, shows the factors influencing gross primary productivity during different construction periods. (a) represents the entire study period, (b) the construction period, (c) combines the pre-construction and construction periods, and (d) combines the construction and post-construction periods. VEG represents vegetation-related variables; Temp represents temperature-related factors; Pre represents precipitation-related factors; and LSI represents landscape pattern-related index. + and - represent positive and negative contributions, respectively.

[0224] Analysis of factors influencing the spatial distribution of GPP revealed that VEG was the primary driver of vegetation GPP in the affected area throughout the study period, accounting for 35.21% of the affected area. This was followed by Temp and Pre, accounting for 24.66% and 25.01%, respectively. LSI had the lowest contribution, at 15.12%. During the construction period, VEG was the dominant factor influencing vegetation GPP in the affected area, accounting for 30.87% of the affected area; this was followed by LSI and Pre, accounting for 27.91% and 21.15%, respectively. Temp dominated the smallest area, at 20.07%. Combining the pre-construction and construction periods, VEG was the primary driver of vegetation GPP in the affected area, accounting for 31.02% of the affected area, followed by LSI and Temp, accounting for 29.69% and 20.02%, respectively. Pre dominated the smallest area, at 19.27%. Combining the construction period and post-construction period, VEG continues to be the dominant factor, accounting for 35.37% of the affected area, followed by Temp and Pre, accounting for 27.04% and 21.83% respectively. LSI dominates the least area, accounting for 15.76%.

[0225] Regarding the changes in the areas of vegetation GPP degradation and improvement, throughout the study period, the areas of vegetation GPP degradation and improvement due to VEG accounted for 29.72% and 5.49% of the affected area, respectively. The areas of vegetation GPP degradation and improvement due to LSI accounted for 6.88% and 8.24% of the affected area, respectively. During the construction period, the areas of degradation and improvement due to VEG accounted for 27.38% and 3.49%, respectively; the areas of degradation and improvement due to LSI accounted for 12.91% and 15%, respectively. Combined before and during construction, the areas of degradation and improvement due to VEG accounted for 26.97% and 4.05%, respectively; the areas of degradation and improvement due to LSI accounted for 14.45% and 15.24%, respectively. Combined during the construction period and after construction, the areas of degradation and improvement due to VEG accounted for 28.96% and 6.41%, respectively, while the areas of degradation and improvement due to LSI accounted for 7.62% and 8.74%, respectively. See Table 5 for details.

[0226] Table 5: Area proportion of each influencing factor in different construction periods

[0227]

[0228] The contribution of each influencing factor to vegetation GPP calculated by the contribution estimation method is as follows: during the entire study period, the contribution of VEG to vegetation GPP is 0.354, the contribution of LSI to vegetation GPP is 0.159, and the contributions of Temp and Pre to vegetation GPP are 0.243 and 0.244 respectively; during the construction period, the contribution of VEG to vegetation GPP is 0.325, the contribution of LSI to vegetation GPP is 0.295, and the contributions of Temp and Pre to vegetation GPP are 0.18 respectively. 7 and 0.193; combined with the pre-construction and construction period, the contribution of VEG to vegetation GPP was 0.306, the contribution of LSI to vegetation GPP was 0.332, and the contributions of Temp and Pre to vegetation GPP were 0.187 and 0.175, respectively; combined with the construction period and post-construction, the contribution of VEG to vegetation GPP was 0.324, the contribution of LSI to vegetation GPP was 0.147, and the contributions of Temp and Pre to vegetation GPP were 0.249 and 0.198, respectively, see Table 6.

[0229] Table 6: Contribution of various influencing factors to vegetation GPP at different construction periods

[0230]

[0231] In other words, VEG is the primary factor controlling GPP across all periods, while the impact of human activities, as measured by LSI, varies with the construction period. The combined LSI before and during construction has the highest degree of control over GPP, indicating that hydropower station construction has altered the basin's ecological environment.

[0232] Figure 12 This is a flow chart of the method for monitoring gross primary productivity of vegetation based on downscaling and multi-source remote sensing data of the present invention, which summarizes the above-mentioned specific embodiments of the present invention. For specific steps and details, please refer to the above description.

[0233] In order to better illustrate the present invention, numerous specific details are provided in the above detailed description. It should be understood by those skilled in the art that the present invention can be practiced without certain specific details. In some instances, methods, means, and components well known to those skilled in the art are not described in detail in order to highlight the main purpose of the present invention.

[0234] The above description is only a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any technician familiar with this technical field can easily think of changes, replacements or omissions within the technical scope disclosed by the present invention, which should be covered by the scope of protection of the present invention.

Claims

1. A method for monitoring gross primary productivity of vegetation based on downscaling and multi-source remote sensing data, characterized in that: The steps include: Step S1, obtaining multi-source remote sensing data of a first resolution, and downscaling the multi-source remote sensing data using a CNN; wherein, in the process of downscaling the multi-source remote sensing data using the CNN, the input layer of the CNN includes at least DEM data and LUCC data, and the background area weight is dynamically adjusted in the loss function using a binary mask of the LUCC data, and in the binary mask of the LUCC data, the vegetation area is 1 and the non-vegetation area is 0; Step S2, simulating vegetation GPP using the CASA method to generate a multi-year GPP dataset for the target area, wherein the maximum light energy utilization rate parameter in the CASA method is determined by the vegetation type; Step S3: For the GPP dataset of the target area over many years, the temporal and spatial variation characteristics of GPP are quantitatively analyzed using the Theil-Sen slope estimation method, the Mann-Kendall non-parametric test method, and the coefficient of variation analysis method, and the spatial distribution and temporal variation pattern of GPP are extracted pixel by pixel; In step S4, principal components of similar factors in environmental variables are extracted according to the principal component analysis method. Then, the local correlation coefficients of each principal component and the GPP are calculated pixel by pixel using the partial correlation method within the GPP range. Finally, the dominant influencing factors of each pixel are identified through the contribution decomposition model, and the environmental factor driving map at the engineering scale is extracted.

2. The method for monitoring gross primary productivity of vegetation based on downscaling and multi-source remote sensing data according to claim 1, characterized in that: The following steps are also included: Step S0: extract the digital river network of the watershed, divide the watershed area, and determine the watershed area controlled by the ridge line as the target area affected by the engineering construction.

3. The method for monitoring gross primary productivity of vegetation based on downscaling and multi-source remote sensing data according to claim 2, characterized in that: The following steps are also included: Using the hydrological analysis module of ArcGIS and the slope runoff simulation method, different regional thresholds were set in the process of extracting the digital river network of the basin.

4. The method for monitoring gross primary productivity of vegetation based on downscaling and multi-source remote sensing data according to claim 3, characterized in that: The extraction of the digital river network of the watershed and the division of the watershed area specifically include: In ArcGIS, the Depression Fill tool of the Spatial Analyst module is called to determine the depression depth through iterative calculation and set the filling threshold. The depression elevation is raised to the lowest value of the adjacent unit to form a depression-free DEM. The flow direction is determined by calculating the distance weight difference between the central grid and the neighboring cells, and the output is a raster data containing 1~255 flow direction codes; Based on the no-sag DEM and flow direction data, the flow accumulation tool is used to count the number of upstream water catchment units of each grid to generate the flow accumulation raster; A catchment area threshold is set to distinguish between river channels and non-river channels, and grids with accumulation greater than or equal to the threshold are identified as river networks. Taking the valley as the pour point, the upstream watershed area is traced back in combination with the flow direction data, the Capture Pour Point tool is called to locate the outlet, and then the Watershed Tool is used to extract all the watershed grids in the basin to generate the watershed boundary extracted by vector.

5. The method for monitoring gross primary productivity of vegetation based on downscaling and multi-source remote sensing data according to claim 4, characterized in that: The following steps are also included: The multi-source remote sensing data is unified to a first resolution by a bilinear interpolation method to obtain multi-source remote sensing data of the first resolution.

6. The method for monitoring gross primary productivity of vegetation based on downscaling and multi-source remote sensing data according to claim 5, characterized in that: The following steps are also included: In the contribution decomposition model, firstly, the spatial distribution data of interannual variation of GPP and environmental factors were obtained through PCA and multi-factor partial correlation analysis; Then, the contribution value estimation method is used to extract the contribution of the interannual variation of environmental factors to the interannual variation of GPP in the target area affected by engineering construction on a pixel-by-pixel basis.

7. The method for monitoring gross primary productivity of vegetation based on downscaling and multi-source remote sensing data according to claim 6, characterized in that: From the patch level, the maximum patch index is selected to represent the spatial proportion of dominant landscape types, and the patch density is selected to reflect the degree of landscape fragmentation. At the type level, the landscape shape index was selected to measure the complexity of patch morphology; At the landscape level, the aggregation index was selected to evaluate the spatial connectivity of patches, and the Shannon diversity index and Shannon evenness index were used to extract the diversity characteristics of the ecosystem.

8. The method for monitoring gross primary productivity of vegetation based on downscaling and multi-source remote sensing data according to claim 7, characterized in that: Also includes: The multi-source remote sensing data includes terrain data and environmental factor data; The environmental variables are divided into four categories: temperature climate variables, moisture-related variables, vegetation characteristic variables and human activity representation variables.

9. The method for monitoring gross primary productivity of vegetation based on downscaling and multi-source remote sensing data according to claim 8, characterized in that: The following steps are also included: The vegetation gross primary productivity monitoring method based on downscaling and multi-source remote sensing data is applied to an ecological remote sensing data processing system.

Images