Global protected area ecological performance quantification assessment method based on time-series remote sensing core variables

Abstract

The application discloses a global protected area ecological effect quantitative evaluation method based on time sequence remote sensing core variables, and belongs to the technical field of ecological protection and remote sensing. The method comprises the following steps: obtaining a target protected area boundary vector, a base period and a final state year multi-source time sequence remote sensing data, a human activity interference data set and a geographic space data set, dividing the protected area and a non-target protected area region in the same region into evaluation grids to respectively form a protected area grid set and a control grid set. System state variables and interference state variables of the two types of grids are obtained through inversion, and corresponding data sets are constructed. Index values of the protection object dimension and the protection effect dimension are calculated based on the data sets, and environmental background interference is removed through the propensity score matching method. After normalization, the weighted analytic hierarchy process is used to obtain a comprehensive score of the protection effect, an effect identifier is generated by comparing with preset grade thresholds, and an evaluation result is output. The application can be adapted to different types of global protected areas, and has high evaluation efficiency and objective and reliable results.

Description

Technical Field

[0001] This invention relates to the fields of ecological protection and remote sensing technology, specifically to a method for quantitatively assessing the ecological effectiveness of global protected areas based on time-series remote sensing core variables. Background Technology

[0002] Current global assessments of the ecological effectiveness of protected areas are mostly based on single-period or low-frequency medium- and low-resolution remote sensing data, with data sampling intervals generally ranging from 5 to 10 years. This makes it difficult to capture interannual ecological dynamics such as shrub invasion, seasonal grazing disturbance, and short-term disaster recovery. Moreover, most assessments focus only on single-dimensional indicators such as the proportion of ecosystem area and the population size of key species, and have not yet established a quantitative assessment system covering the three levels of ecosystem, species, and genes. In particular, they lack the inclusion of core indicators that represent long-term evolutionary heritage, such as phylogenetic diversity, and cannot fully reflect the effectiveness of protected areas in maintaining the entire spectrum of biodiversity.

[0003] Current assessments largely rely on manual field sampling surveys, which are highly subjective in sample selection and have limited coverage. Different countries and protected areas employ significantly different indicator calculation methods and data sources. Some institutions use simple weighting methods, while others use non-standardized expert scoring methods, resulting in inconsistent normalization rules. This leads to a lack of comparability of assessment results across regions and scales, making it difficult to support horizontal benchmarking of global protected area effectiveness. Traditional assessments often use a single-sample analysis paradigm of the base-to-final state difference within the protected area, failing to include control samples from non-protected areas in the same location. They also fail to match key covariates affecting ecological conditions, such as altitude, slope, distance from roads and settlements, and cannot isolate background interferences such as climate change, interannual precipitation fluctuations, and development activities in surrounding areas. For example, the increase in vegetation cover in a subtropical protected area in Southwest China may simply be due to abnormally high precipitation during the assessment period, rather than the effectiveness of management measures. Current technology cannot distinguish between the net increase brought about by protection measures and the spurious effects of changes in the natural background. In addition, existing assessments rely heavily on localized data preprocessing processes and do not integrate publicly shared resources such as Sentinel, Landsat time-series remote sensing images, and the Global Human Footprint dataset of human activity disturbances. This results in high data acquisition costs, low processing efficiency, and a single assessment cycle that is often 3-5 years behind, making it impossible to support dynamic monitoring and adaptive management decisions for protected areas.

[0004] Therefore, there is an urgent need for a quantitative assessment method for the ecological effectiveness of global protected areas that is based on multi-source time-series remote sensing big data, can eliminate environmental background interference, and has a high degree of standardization. Summary of the Invention

[0005] The purpose of this invention is to provide a quantitative assessment method for the ecological effectiveness of global protected areas based on time-series remote sensing core variables, in order to solve the problems of existing assessment methods being unable to eliminate environmental background interference, having low standardization, and being difficult to compare across regions.

[0006] To achieve the above objectives, the technical solution adopted by the present invention is as follows: A quantitative assessment method for the ecological effectiveness of global protected areas based on time-series remote sensing core variables includes the following steps: Acquire boundary vector data of the target protected area and multi-source time-series remote sensing satellite imagery data of the base year and the final state year, as well as publicly available human activity disturbance datasets and geospatial datasets; The spatial range of the target protected area is defined based on the boundary vector data, and the target protected area and the corresponding area of ​​the non-target protected area in the same location are divided into several evaluation grids according to a preset spatial resolution. The evaluation grids located within the target protected area constitute the protected area grid set, and the evaluation grids located in the corresponding area of ​​the non-target protected area in the same location constitute the control grid set. By combining the aforementioned time-series remote sensing satellite image data, vegetation indices, surface temperature, and leaf area indices of each grid in the protected area grid set and control grid set in the base year and the final state year are obtained as ecosystem state variables. Land cover change and traces of human activities based on remote sensing image interpretation are used as disturbance state variables. A gridded assessment dataset containing remote sensing spectral characteristics, land cover type, environmental covariates, and species distribution information is constructed. Data is then extracted from the dataset to construct the protected area grid set dataset and the control grid set dataset, respectively. Based on the protected area grid dataset, the index values ​​of the protected object dimension are calculated with the target protected area as a whole. The index values ​​of the protected object dimension include the natural ecosystem integrity index, the protected coverage rate of threatened species habitat, and the phylogenetic diversity protection index, which correspond to the evaluation indicators at the ecosystem level, species level, and gene level, respectively. Based on the protected area grid dataset and the control grid dataset, the index values ​​of the protection effectiveness dimension are calculated. The index values ​​of the protection effectiveness dimension include the net change rate of natural vegetation cover, the net change rate of key species populations, and the net reduction rate of human activity disturbance, which correspond to the assessment indicators at the ecosystem level, species level, and threat mitigation level, respectively. The net change rate of natural vegetation cover and the net reduction rate of human activity disturbance are matched with each grid in the control grid dataset using the propensity score matching method, and the difference in the index change rate between grids is used as the net change after removing environmental background disturbance. The index values ​​of the protected object dimension and the index values ​​of the protection effectiveness dimension are normalized, and then weighted and summed according to the preset weights determined by the analytic hierarchy process based on the consistency test of the judgment matrix to obtain the comprehensive score of the protection effectiveness of the target protected area. Based on the comprehensive score of the protection effectiveness and the preset level threshold table, the protection effectiveness level identifier is automatically generated, and the quantitative assessment result of the ecological effectiveness of the target protected area is output.

[0007] Furthermore, the area corresponding to the non-target protected area in the same location is the area covered by a preset buffer zone extending outward from the boundary of the target protected area.

[0008] Furthermore, the index values ​​for the protected object dimension include the natural ecosystem integrity index, the habitat protection coverage rate of threatened species, and the phylogenetic diversity protection index, specifically: The natural ecosystem integrity index is used to assess the area proportion and spatial pattern integrity of natural ecosystems within the target protected area. It is calculated based on natural ecosystem patch data in the protected area grid dataset, which is obtained by supervised classification of the multi-source time-series remote sensing satellite image data. The coverage rate of threatened species habitat protection is used to assess the coverage of the target protected area on the suitable habitat of threatened species. It is based on the vegetation type, altitude and climate variable data in the protected area grid dataset to generate a suitable habitat distribution map, and then the suitable habitat distribution map is spatially overlaid with the boundary of the target protected area and the average value is obtained by weighting according to the species endangerment level. The phylogenetic diversity conservation index is used to assess the proportion of the total evolutionary history represented by the species set within the target protected area relative to the overall regional population. It is calculated based on a pre-constructed regional species phylogenetic tree and species distribution data in the protected area grid set dataset. The region includes the target protected area and the corresponding region of the non-target protected area in the same location. The phylogenetic diversity conservation index is the ratio of the sum of the phylogenetic diversity values ​​of each grid in the protected area grid set within the region to the sum of the phylogenetic diversity values ​​of each grid in the protected area grid set and the control grid set within the region.

[0009] Furthermore, the indicator values ​​for the conservation effectiveness dimension include the net change rate of natural vegetation cover, the net change rate of key species populations, and the net reduction rate of human activity disturbance, specifically: The calculation process of the net change rate of natural vegetation cover is as follows: extract the natural vegetation area of ​​the protected area grid set dataset and the control grid set dataset in the base year and the final year, respectively. The natural vegetation area is obtained by statistically analyzing the area of ​​the forest, grassland and wetland types in the natural ecosystem patch data. Calculate the area change rate of each grid. Subtract the area change rate of each grid in the control grid set from the area change rate of each grid in the protected area grid set. The difference is the net change rate. The calculation process for the net population change rate of the key species is as follows: Population monitoring data of the key protected species in the base year and terminal year are obtained from the protected area grid dataset and the control grid dataset. The population monitoring data includes at least one of the relative abundance index monitored by infrared cameras, transect survey encounter rate, or fixed quadrat count. The population change rate of each grid in the protected area grid dataset and each grid in the control grid dataset is calculated from the base year to the terminal year. The difference between the population change rate of each grid in the protected area grid dataset and the population change rate of each grid in the control grid dataset is the net population change rate of the key species. The calculation process of the net reduction rate of human activity interference is as follows: based on the interpretation of the multi-source time-series remote sensing satellite image data, or based on the human activity interference intensity obtained from the human activity interference dataset, the human activity interference intensity of the protected area grid set dataset and the control grid set dataset in the base year and the final year are obtained respectively, the rate of change of interference intensity of each is calculated, and then the rate of change of each grid in the control grid set is subtracted from the rate of change of each grid in the protected area grid set, and the difference is the net reduction rate.

[0010] Furthermore, the propensity score matching method specifically includes: Based on the geospatial dataset, environmental variables of each grid in the protected area grid set and the control grid set are extracted simultaneously as matching covariates. The environmental variables include at least one of altitude, slope, aspect, distance from road, distance from settlement, annual precipitation and average annual temperature. Using a logistic regression algorithm, with the matching covariate as the independent variable and the group to which the evaluation grid belongs as the dependent variable, the tendency score of each grid in the protected area grid set and the control grid set to belong to the protected area grid set is calculated. The nearest neighbor matching method is used to identify the grid in the control grid set that has the smallest absolute difference in matching tendency scores among all grids in the protected area grid set. The matching results are tested for balance. The standardized mean differences of each matching covariate in the matched protected area grid set and the control grid set are calculated. When the absolute value of all standardized mean differences is less than 0.25, the matching results are deemed to meet the covariate distribution consistency requirements.

[0011] Furthermore, the normalization process employs the Min-Max normalization method, linearly mapping all original index values ​​to the [0,1] interval. Positive indices use a forward normalization formula, while negative indices use a reverse normalization formula, ensuring that all indices are positive after standardization. The larger the value, the better the protection effect.

[0012] Furthermore, the construction process of the preset weights determined by the analytic hierarchy process based on the consistency test of the judgment matrix specifically includes: The weights are general weights that are pre-constructed based on knowledge in the field of ecological protection and can be reused for different target protected areas; A three-layer hierarchical structure is constructed, comprising a target layer, a criterion layer, and an indicator layer. The target layer is the comprehensive score of the protection effectiveness, the criterion layer includes the protection object dimension and the protection effectiveness dimension, and the indicator layer includes all indicator values ​​of the protection object dimension and the protection effectiveness dimension. Based on the preset 1-9 scale method, a criterion layer judgment matrix and an indicator layer judgment matrix are constructed. The elements of the judgment matrix represent the relative importance values ​​of any two indicators within the same level. Calculate the maximum eigenvalue and corresponding eigenvector of each judgment matrix, and normalize the eigenvectors to obtain the weight vectors of each dimension of the criterion layer and the weight vectors of each index value of the index layer. Calculate the consistency ratio of each judgment matrix. When the consistency ratio is less than 0.1, the weight vector is deemed valid.

[0013] Furthermore, the preset level threshold table includes four levels, specifically: Excellent level, corresponding to a comprehensive score of no less than 85 points for the protection effectiveness; A good rating corresponds to a comprehensive score of protection effectiveness of not less than 70 points and less than 85 points. Medium level, corresponding to a comprehensive score of protection effectiveness of not less than 60 points and less than 70 points; A poor rating corresponds to a comprehensive protection effectiveness score of less than 60 points.

[0014] Furthermore, the results of the quantitative assessment of ecological effectiveness are presented through at least one of the following visualization methods: Radar charts are used to display the standardized score distribution of each indicator value in the dimensions of protected objects and protection effectiveness of the target protected area, and to identify shortcomings in protection effectiveness. Spatial heat map is used to map the protection effectiveness dimension index values ​​of each grid in the grid set of the target protected area to the spatial range of the target protected area, and to identify the spatial heterogeneity of protection effectiveness; The traffic light warning map is used to mark the changing trends of each indicator value of the target protected area in terms of the protection object dimension and the protection effectiveness dimension from the base period to the final state year in three categories: red light indicates that the indicator value has decreased compared with the base year, yellow light indicates that the change in the indicator value compared with the base year is within a preset threshold, and green light indicates that the indicator value has increased compared with the base year.

[0015] Furthermore, the multi-source time-series remote sensing satellite image data originates from at least one of the Sentinel series, Landsat series, or MODIS remote sensing satellite products; The human activity disturbance dataset includes the Global Human Footprint dataset; The geospatial dataset includes at least one of digital elevation models, road network data, settlement distribution data, and climate datasets.

[0016] Compared with the prior art, the present invention has the following beneficial effects: 1. This invention adopts a two-dimensional and three-level panoramic assessment framework, which organically combines the three biological levels of ecosystem, species and genes with the two management dimensions of the value of the protected object and the effectiveness of protection management, thus solving the problem of quantifying the effectiveness of evolutionary diversity protection in global protected area assessment.

[0017] 2. Introducing counterfactual causal inference to eliminate environmental noise. The propensity score matching method is applied to the assessment of protected area effectiveness. By finding control grids with similar environmental backgrounds in non-protected areas within the same location for each protected area grid, the difference in the rate of change of indicators is calculated, effectively eliminating environmental background interference such as climate change and macroeconomic fluctuations, ensuring that the assessment results truly reflect the net effect of protection measures.

[0018] 3. High degree of standardization and automation. Based on cloud computing platforms such as Google Earth Engine (GEE), and integrating public datasets such as Sentinel, Landsat, MODIS, and Global Human Footprint, it has achieved full-process automation from data acquisition, grid computing, indicator fusion to grade output, and supports annual dynamic assessment of protected areas at the global scale.

[0019] 4. By using Min-Max normalization and AHP (Analytical Hierarchy Process), multidimensional indicators are integrated into a single comprehensive score, which is then visualized through radar charts, spatial heat maps, and traffic light warning maps. This allows for both macroscopic control of the overall condition of protected areas and precise location of ecological degradation and weak indicators, providing a scientific basis for adaptive management. Attached Figure Description

[0020] Figure 1 This is a flowchart of the method of the present invention. Detailed Implementation

[0021] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention. Furthermore, the technical features involved in the various embodiments of this invention described below can be combined with each other as long as they do not conflict with each other.

[0022] This embodiment provides a method for quantitatively assessing the ecological effectiveness of global protected areas based on time-series remote sensing core variables. The main process is as follows: Figure 1 As shown, this implementation utilizes the Google Earth Engine (GEE) cloud computing platform. The assessment object in this embodiment is a mid-subtropical mountain forest nature reserve, primarily protecting the mid-subtropical primary forest ecosystem and threatened wild flora and fauna distributed within the region. The assessment period is set at 10 years according to regulations, a time window consistent with the standard periodic requirements of the UN Kunming-Montreal Global Biodiversity Framework for global protected area effectiveness assessments from 2011 to 2020. Base year. For 2010, the final state year. For 2020, both data sets used publicly available Landsat / Sentinel time-series remote sensing products, eliminating the need for additional monitoring time.

[0023] Step 1: Acquire the boundary vector data of the target protected area and multi-source time-series remote sensing satellite imagery data for the base year and the final state year, as well as publicly available human activity disturbance datasets and geospatial datasets. The multi-source time-series remote sensing satellite imagery data comes from at least one of the Sentinel series, Landsat series, or MODIS remote sensing satellite products; the human activity disturbance dataset includes the Global Human Footprint dataset; the geospatial dataset includes at least one of digital elevation models, road network data, settlement distribution data, and climate datasets.

[0024] In this embodiment, the multi-source time-series remote sensing satellite imagery data comes from Sentinel-2 MSI (10m resolution) and Landsat 7 ETM+ (30m resolution) remote sensing products; the publicly available human activity disturbance dataset uses the Global Human Footprint (GHF) dataset; the geospatial dataset includes the SRTM digital elevation model, OpenStreetMap road network data, WorldPop settlement distribution data, and WorldClim v2.1 climate dataset.

[0025] Step 2: Define the spatial range of the target protected area based on the boundary vector data, and divide the target protected area and the corresponding area of ​​the non-target protected area in the same location into several evaluation grids according to a preset spatial resolution. The evaluation grids located within the target protected area constitute the protected area grid set, and the evaluation grids located in the corresponding area of ​​the non-target protected area in the same location constitute the control grid set.

[0026] Specifically, the target protected area and the corresponding non-target protected area in the same location are divided into 302 assessment grids with a preset spatial resolution of 1km × 1km. According to industry standards, a 1km spatial resolution is the standard grid scale for global protected area ecological assessment, effectively characterizing the spatial heterogeneity of the ecosystem while maintaining computational efficiency. Assessment grids located within the target protected area constitute the protected area grid set, while assessment grids located within the corresponding non-target protected area in the same location constitute the control grid set. The corresponding non-target protected area in the same location is the area covered by a preset buffer zone extending outward from the boundary of the target protected area, such as 10km.

[0027] Step 3: Combine the time-series remote sensing satellite image data to retrieve the vegetation index, surface temperature, and leaf area index of each grid in the protected area grid set and the control grid set in the base year and the final state year as ecosystem state variables, and land cover change and traces of human activities based on remote sensing image interpretation as disturbance state variables. Construct a gridded assessment dataset containing remote sensing spectral characteristics, land cover type, environmental covariates, and species distribution information, and extract data from it to construct the protected area grid set dataset and the control grid set dataset respectively.

[0028] In one specific embodiment, ecosystem state variables include Normalized Difference Vegetation Index (NDVI), Leaf Area Index (LAI), and Land Surface Temperature (LST). NDVI / EVI is retrieved based on the Sentinel-2 red-edge band, LAI is retrieved based on the PROSAIL radiative transfer model, and LST is retrieved based on the thermal infrared band. Disturbance state variables include land cover change information detected using the CCDC algorithm and traces of human activities, such as construction land expansion and farmland encroachment, based on GHF datasets and remote sensing interpretation. Finally, the aforementioned gridded assessment dataset is constructed, and data is extracted from it to construct protected area grid datasets and control grid datasets, respectively.

[0029] Step 4: Based on the protected area grid dataset, calculate the indicator values ​​for the protected object dimension, taking the target protected area as a whole. The indicator values ​​for the protected object dimension include the natural ecosystem integrity index, the habitat protection coverage rate for threatened species, and the phylogenetic diversity protection index, corresponding to assessment indicators at the ecosystem level, species level, and genetic level, respectively. These include: (1) The Natural Ecosystem Integrity Index, i.e., the ecosystem hierarchy, is used to assess the area proportion and spatial pattern integrity of natural ecosystems within the target protected area. It is calculated based on natural ecosystem patch data in the protected area grid dataset, which is obtained through supervised classification of the multi-source time-series remote sensing satellite image data, covering three types of natural ecosystems: forest, grassland, and wetland. In this embodiment, a random forest classifier is used to implement supervised classification. The training samples are from the publicly available Global Land Cover SHARE (GLC_SHARE) dataset, and the overall classification accuracy is better than 85%, which meets the conventional requirements for ecological remote sensing classification.

[0030] First, calculate the area ratio of natural ecosystems: in, The percentage of the area covered by natural ecosystems (%) Let be the area (km²) of the j-th natural ecosystem patch; The total number of patches in the natural ecosystem; The total area of ​​the protected area (km²).

[0031] Calculate the fragmentation index of natural ecosystems: in, It is the natural ecosystem fragmentation index, with a value range of [0, 1]. The larger the value, the more severe the fragmentation.

[0032] Calculate the area-weighted fractal dimension: in, is the area-weighted fractal dimension, with values ​​ranging from [1, 2]. Let be the perimeter (km) of the j-th natural ecosystem patch, obtained by calculating the geometric boundary length after raster-to-vector conversion; lg is a common logarithm with base 10.

[0033] The final comprehensive index is the natural ecosystem integrity index: in, The value ranges from [0, 100], with larger values ​​indicating a more complete spatial pattern of the ecosystem. This embodiment calculates... =81.

[0034] (2) Habitat protection coverage of threatened species, i.e., species level, is used to assess the coverage of suitable habitats for threatened species by the target protected area. Based on vegetation type, altitude, and climate variables in the protected area grid dataset, a distribution map of suitable habitats for threatened species in the region is generated using the MaxEnt model. The input environmental variables of the MaxEnt model include, but are not limited to, Normalized Difference Vegetation Index (NDVI), Leaf Area Index (LAI), Land Surface Temperature (LST), altitude, slope, and annual precipitation. The distribution map is overlaid with the protected area boundary for analysis, and the average coverage is calculated by weighting according to the species endangerment level (CR level weight 3, EN level weight 2, VU level weight 1). The threatened species include VU vulnerable, EN endangered, and CR critically endangered, and are sourced from the IUCN Red List and the National Key Protected Wild Animals and Plants List.

[0035] Calculate the coverage of a single species: Among them, in: Habitat protection coverage (%) for species i; , where is the area (km²) of the suitable habitat of species i that is covered by protected areas. The total suitable habitat area (km²) of species i within the study area (including protected areas and surrounding buffer zones).

[0036] Calculate the overall coverage rate: in, The overall coverage rate of habitat protection for threatened species (%). is the endangerment weight of the i-th species, which can be assigned different weights according to the IUCN endangerment level, such as CR=3, EN=2, VU=1; m is the total number of threatened species.

[0037] The calculation in this embodiment is =87%.

[0038] (3) The phylogenetic diversity conservation index, i.e., at the gene level, is used to assess the proportion of the total evolutionary history represented by the species set within the target protected area relative to the overall regional population. It is calculated based on species distribution data from a pre-constructed regional species phylogenetic tree and a protected area grid dataset. The region includes the target protected area and the corresponding region of the non-target protected area in the same location, i.e., the union of the target protected area boundary with a 10km pre-defined buffer zone. The regional species phylogenetic tree covers vascular plant species within the protected area and its surrounding pre-defined buffer zone, and is constructed using the Open Tree of Life database. The calculation process for the phylogenetic diversity conservation index is as follows: Developmental diversity values ​​for single-grid phylogenetic systems: in, Let k be the phylogenetic diversity value of grid k; This refers to all branches in the phylogenetic tree involving the species set within grid k; The evolutionary time length of branch b (millions of years) is derived from the branch length calibration results of the Open Tree of Life database.

[0039] Total PD value of protected areas: PA represents all grids within the protected area grid set.

[0040] Total Regional PD: in, This includes all grids within the union of the protected area and the preset buffer zone, including the protected area grid set and the control grid set.

[0041] Phylogenetic diversity protection index: in, The phylogenetic diversity conservation index (%) reflects the proportion of evolutionary history captured in protected areas relative to the total evolutionary history of the region. In this example, the IPD was calculated to be 67%.

[0042] Step 5: Based on the protected area grid dataset and the control grid dataset, calculate the indicator values ​​for the protection effectiveness dimension. These indicator values ​​include the net change rate of natural vegetation cover, the net change rate of key species populations, and the net reduction rate of human activity disturbance, corresponding to assessment indicators at the ecosystem, species, and threat mitigation levels, respectively. The net change rate of natural vegetation cover and the net reduction rate of human activity disturbance are matched to each grid in the protected area grid dataset using a propensity score matching method, with the difference in indicator change rates between grids used as the net change after removing environmental background disturbances. The specific process is as follows: (1) Net change rate of natural vegetation cover, i.e., ecosystem level. The natural vegetation area of ​​the protected area grid set dataset and the control grid set dataset in the base year and the final year were extracted respectively. The natural vegetation area was obtained by statistically analyzing the area of ​​the forest, grassland and wetland types in the natural ecosystem patch data. The area change rate of each was calculated. The net change rate was obtained by subtracting the area change rate of each grid in the control grid set from the area change rate of each grid in the protected area grid set.

[0043] Rate of change of protected areas: , The areas of natural vegetation (km²) within the protected area grid set in the base period and final state years are respectively.

[0044] Comparison grid change rate: in, , The figures represent the natural vegetation area (km²) within the baseline and final state control grid sets, respectively.

[0045] Net change rate: NRV represents the net change rate of natural vegetation cover (%). A positive value indicates that the vegetation restoration effect in the protected area is better than that in the control area, while a negative value indicates that the vegetation restoration effect in the protected area is worse than that in the external control area. In this example, NRV was calculated to be 2.8%.

[0046] (2) Net population change rate of key species, i.e., species level. Population monitoring data of key protected species in the base year and terminal year were obtained from the protected area grid dataset and the control grid dataset. The population monitoring data included at least one of the Relative Abundance Index (RAI) from infrared camera monitoring, transect survey encounter rate, or fixed quadrat count. The population change rate of each grid in the protected area grid dataset and each grid in the control grid dataset was calculated from the base year to the terminal year. The difference between the population change rate of each grid in the protected area grid dataset and the population change rate of each grid in the control grid dataset was the net population change rate of the key species. The Relative Abundance Index (RAI) was calculated as: number of independent valid photos of a species taken by an infrared camera / (camera working days × number of cameras), which is a wildlife population monitoring algorithm.

[0047] Population change rate in protected areas: in, , These represent the population size or relative abundance index of key species within the protected area grid set during the base and final states, respectively. Population change rate in the control area: in, , These represent the population size or relative abundance index of key species within the control grid sets for the base and final states, respectively. Net change rate of a single species: in, Let be the net change rate (%) of species i, where i is the i-th key protected species; Overall net change rate: Where k represents the total number of key protected species. This embodiment calculates... =7.2%.

[0048] (3) Net reduction rate of human activity disturbance, i.e. threat mitigation level. Based on the interpretation of multi-source time-series remote sensing satellite image data, or based on the human activity disturbance intensity obtained from the human activity disturbance dataset GHF, the human activity disturbance intensity of the protected area grid set dataset and the control grid set dataset in the base year and the final year are obtained respectively. The rate of change of disturbance intensity of each is calculated. Then, the rate of change of each grid in the control grid set is subtracted from the rate of change of each grid in the protected area grid set. The difference is the net reduction rate.

[0049] Rate of change of disturbance in protected areas: in, , These represent the intensity of human activity disturbance within the protected area grid set during the base and final states, respectively. Comparison with grid disturbance change rate: in, , These represent the intensity of human activity disturbance within the control grid sets for the base and final states, respectively. Net reduction rate: Wherein, NHR represents the net reduction rate of human disturbance (%). A positive value indicates that the rate of disturbance growth in the protected area is slower than that in the control area, indicating significant protection effectiveness; the higher the value, the better the effectiveness. A negative value indicates that the rate of disturbance growth in the protected area is slower than that in the control area. In this embodiment, NHR was calculated to be 11.6%.

[0050] The control grid sets for the net change rate of natural vegetation cover and the net reduction rate of human activity disturbance mentioned above were constructed using propensity score matching, specifically as follows: Based on the geospatial dataset, environmental variables of each grid in the protected area grid set and the control grid set are extracted simultaneously as matching covariates. The environmental variables include at least one of altitude, slope, aspect, distance from road, distance from settlement, annual precipitation and average annual temperature. Using a logistic regression algorithm, with the matching covariate as the independent variable and the group to which the evaluation grid belongs as the dependent variable, the tendency score of each grid in the protected area grid set and the control grid set to belong to the protected area grid set is calculated. The nearest neighbor matching method is adopted, which selects at least one grid from the control grid set with the smallest difference in the absolute value of the matching tendency score among all grids in the protected area grid set. The preferred matching ratio is 1:1 to reduce the estimation variance. If the number of grids in the treatment group is small or the matching difficulty is high, it can be appropriately relaxed to 1:N (N≤5) to improve the matching success rate.

[0051] A balance test was performed on the matching results. The standardized mean differences of each matching covariate in the matched protected area grid set and the control grid set were calculated. When the absolute value of all standardized mean differences was less than 0.25, the matching results were deemed to meet the covariate distribution consistency requirement. According to the general standard of causal inference analysis, a standardized mean difference of <0.25 is the conventional significance threshold for propensity score matching, indicating that there is no systematic bias in the covariate distribution between the treatment group and the control group after matching.

[0052] If the matching results do not meet the consistency requirements, the combination of matching covariates or the matching tolerance can be adjusted, and the matching process can be re-executed until the requirements are met.

[0053] Step 6: Normalize the index values ​​of the protected object dimension and the protected effectiveness dimension, and perform weighted summation according to the preset weights determined by the analytic hierarchy process based on the consistency test of the judgment matrix to obtain the comprehensive score of the protected area's protection effectiveness.

[0054] Normalization was performed using the Min-Max normalization method, applying the three index values ​​of the protection object dimension obtained above ( , , ) and the three indicators of protection effectiveness (NRV, The original values ​​of 6 indicators (NHR) were standardized, and all original values ​​of indicators were linearly mapped to the [0,1] interval.

[0055] All six indicators involved in this embodiment are positive indicators; the larger the value, the better the protection effect. A positive normalization formula is used to ensure that all indicators are positive after standardization. in, and These are the minimum and maximum values ​​of the corresponding indicators in all evaluation grids within the study area, including both the target protected area and the corresponding non-target protected area in the same location, thus reflecting the relative level of the indicators in the regional context.

[0056] If inverse indicators are involved (such as the original value of human activity disturbance intensity), then the inverse normalization formula is used: Subsequently, the weights were determined using a pre-defined analytic hierarchy process based on the consistency test of the judgment matrix. These weights are general weights pre-constructed based on knowledge in the field of ecological protection and can be reused for different target protected areas. The construction process is as follows: A three-tiered hierarchical structure is constructed, comprising an objective layer, a criterion layer, and an indicator layer. The objective layer represents the comprehensive score of protection effectiveness; the criterion layer includes the dimensions of protected objects and protection effectiveness; and the indicator layer includes three specific indicators for both the protected object dimension and the protection effectiveness dimension. Based on the pre-defined 1-9 scale method, and combined with prior knowledge in the field of ecological assessment, a judgment matrix for the criterion layer and a judgment matrix for the indicator layer are constructed. The elements of the judgment matrix represent the relative importance values ​​of any two indicators within the same level. The maximum eigenvalue and corresponding eigenvector of each judgment matrix are calculated, and the eigenvectors are normalized to obtain the weight vectors of each dimension of the criterion layer and the weight vectors of each index value of the index layer. The typical weight allocation used in this embodiment is as follows: Global weight 0.4 for the protected object dimension, with a weight of 0.5 within the natural ecosystem integrity index dimension and a global weight of 0.2; weight of 0.3 within the threatened species habitat protection coverage rate dimension and a global weight of 0.12; and weight of 0.2 within the phylogenetic diversity protection index dimension and a global weight of 0.08. Global weight 0.6 for the protection effectiveness dimension, with a weight of 0.3 within the net change rate of natural vegetation cover dimension and a global weight of 0.18; weight of 0.4 within the net change rate of key species populations dimension and a global weight of 0.24; and weight of 0.3 within the net reduction rate of human activity disturbance dimension and a global weight of 0.18. The consistency ratio of each judgment matrix is ​​calculated. In this embodiment, the consistency ratio is 0.07, which is less than 0.1, indicating that the weight vector is valid. A consistency ratio of <0.1 is a conventional acceptable threshold for the analytic hierarchy process (AHP), indicating that the logical consistency of the judgment matrix meets the evaluation requirements, and is a commonly used standard for weight validity testing.

[0057] If the consistency ratio is not less than 0.1, the values ​​of the elements in the judgment matrix can be adjusted, and the weight calculation steps can be re-executed until the requirements are met.

[0058] Finally, the standardized scores of the six indicators are weighted and summed according to the above weights to obtain the comprehensive score of the protection effectiveness of the target protected area. The formula is as follows: in, The overall score for protection effectiveness is 100 points. Let be the weight of the j-th criterion layer; Let be the weight of the i-th tertiary indicator in the j-th dimension; Let be the standardized score of the i-th tertiary indicator under the j-th dimension. Substituting the data from this example, we obtain... =76 points.

[0059] Step 7: Based on the comprehensive score of the protection effectiveness and the preset level threshold table, automatically generate the protection effectiveness level label and output the quantitative evaluation result of the ecological effectiveness of the target protected area.

[0060] The preset level threshold table includes four levels, as shown in Table 1.

[0061] Table 1. Preset Level Threshold Table (Classification) In this embodiment, based on the comprehensive protection effectiveness score of 76 points, and referring to the preset level threshold table, this score belongs to the good level, which indicates that the value of the protected object is high and the overall management is effective, and an "good" protection effectiveness level label is automatically generated.

[0062] The results of the quantitative assessment of ecological effectiveness are presented in at least one of the following visualization formats: (1) Radar chart, used to display the standardized score distribution of each indicator value of the protection object dimension and protection effectiveness dimension of the target protected area. In this embodiment, the system development diversity protection index is identified as the current shortcoming indicator of protection effectiveness through radar chart, providing a basis for targeted management and realizing the function of identifying the shortcoming of protection effectiveness. (2) Spatial heat map, used to map the protection effectiveness dimension index value of each grid in the grid set of the target protected area to the spatial range of the target protected area. In this embodiment, the spatial heat map identifies that the net vegetation change rate in the northern area of ​​the protected area is low, which is an ecologically weak area, thus realizing the function of identifying the spatial heterogeneity of protection effectiveness. (3) Traffic light warning map, used to mark the changing trends of each indicator value of the target protected area in the dimensions of protected objects and protection effectiveness from the base period to the end state year in three categories: red light indicates that the indicator value decreases compared with the base year. In this embodiment, the change rate of the phylogenetic diversity index is -2%, and it is marked as red light; yellow light indicates that the change of the indicator value compared with the base year is within the preset threshold. The threshold is set to ±2% with reference to the 2% interannual fluctuation significance test standard commonly used in the field of ecological monitoring; in this embodiment, the change of the net reduction rate of human activity disturbance is 1.2%, which is within the preset threshold of ±2%, and it is marked as yellow light; green light indicates that the indicator value increases compared with the base year. In this embodiment, the net change rate of natural vegetation cover and the net change rate of key species populations both show an upward trend, and it is marked as green light.

[0063] The method described in this embodiment possesses high versatility and scalability, applicable to the assessment of various protected areas across different climate zones and with different protected objects worldwide. Since the core steps are based on standardized multi-source remote sensing products and publicly available geographic datasets, and a unified data structure is achieved through a gridded assessment dataset, this method can overcome the limitations of localized manual survey data and achieve cross-regional reuse. Furthermore, all computational steps can be encapsulated in the GEE cloud computing platform for automated batch processing. Compared to traditional manual ground surveys, this reduces the single assessment cycle from several years to 1-2 months, significantly improving the response speed and timeliness of protected area management, and meeting the timeliness requirements of the "Technical Specifications for the Assessment of the Conservation Effectiveness of Nature Reserves."

[0064] A complete assessment report can be prepared based on the above assessment results, including the assessment background, indicator calculation process, comprehensive score, grade rating, problem analysis and management recommendations, which directly supports the adaptive management of protected areas and the progress assessment of the "3030" goals of the Kunming-Montreal Global Biodiversity Framework.

[0065] Embodiments of the present invention may be provided as methods, systems, or computer program products. Therefore, the present invention may take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention may take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0066] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0067] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0068] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0069] Contents not described in detail in this specification are prior art known to those skilled in the art. It is hereby indicated that the above description is intended to help those skilled in the art understand this invention, but does not limit the scope of protection of this invention. Any equivalent substitutions, modifications, improvements, or simplifications of the above descriptions that do not depart from the essential content of this invention fall within the scope of protection of this invention.

Claims

1. A method for quantitatively assessing the ecological effectiveness of global protected areas based on time-series remote sensing core variables, characterized in that... Includes the following steps: Acquire boundary vector data of the target protected area and multi-source time-series remote sensing satellite imagery data of the base year and the final state year, as well as publicly available human activity disturbance datasets and geospatial datasets; The spatial range of the target protected area is defined based on the boundary vector data, and the target protected area and the corresponding area of ​​the non-target protected area in the same location are divided into several evaluation grids according to a preset spatial resolution. The evaluation grids located within the target protected area constitute the protected area grid set, and the evaluation grids located in the corresponding area of ​​the non-target protected area in the same location constitute the control grid set. By combining the aforementioned time-series remote sensing satellite image data, vegetation indices, surface temperature, and leaf area indices of each grid in the protected area grid set and control grid set in the base year and the final state year are obtained as ecosystem state variables. Land cover change and traces of human activities based on remote sensing image interpretation are used as disturbance state variables. A gridded assessment dataset containing remote sensing spectral characteristics, land cover type, environmental covariates, and species distribution information is constructed. Data is then extracted from the dataset to construct the protected area grid set dataset and the control grid set dataset, respectively. Based on the protected area grid dataset, the index values ​​of the protected object dimension are calculated with the target protected area as a whole. The index values ​​of the protected object dimension include the natural ecosystem integrity index, the protected coverage rate of threatened species habitat, and the phylogenetic diversity protection index, which correspond to the evaluation indicators at the ecosystem level, species level, and gene level, respectively. Based on the protected area grid dataset and the control grid dataset, the index values ​​of the protection effectiveness dimension are calculated. The index values ​​of the protection effectiveness dimension include the net change rate of natural vegetation cover, the net change rate of key species populations, and the net reduction rate of human activity disturbance, which correspond to the assessment indicators at the ecosystem level, species level, and threat mitigation level, respectively. The net change rate of natural vegetation cover and the net reduction rate of human activity disturbance are matched with each grid in the control grid dataset using the propensity score matching method, and the difference in the index change rate between grids is used as the net change after removing environmental background disturbance. The index values ​​of the protected object dimension and the index values ​​of the protected effectiveness dimension are normalized, and then weighted and summed according to the preset weights determined by the analytic hierarchy process based on the consistency test of the judgment matrix to obtain the comprehensive score of the protection effectiveness of the target protected area. Based on the comprehensive score of the protection effectiveness and the preset level threshold table, the protection effectiveness level identifier is automatically generated, and the quantitative assessment result of the ecological effectiveness of the target protected area is output.

2. The method for quantitatively evaluating the ecological effectiveness of global protected areas based on time-series remote sensing core variables as described in claim 1, characterized in that, The area corresponding to the non-target protected area in the same location is the area covered by a preset buffer zone extending outward from the boundary of the target protected area.

3. The method for quantitatively assessing the ecological effectiveness of global protected areas based on time-series remote sensing core variables as described in claim 2, characterized in that, The index values ​​for the protected objects dimension include the natural ecosystem integrity index, the habitat protection coverage rate of threatened species, and the phylogenetic diversity protection index, specifically: The natural ecosystem integrity index is used to assess the area proportion and spatial pattern integrity of natural ecosystems within the target protected area. It is calculated based on natural ecosystem patch data in the protected area grid dataset, which is obtained by supervised classification of the multi-source time-series remote sensing satellite image data. The coverage rate of threatened species habitat protection is used to assess the coverage of the target protected area on the suitable habitat of threatened species. It is based on the vegetation type, altitude and climate variable data in the protected area grid dataset to generate a suitable habitat distribution map, and then the suitable habitat distribution map is spatially overlaid with the boundary of the target protected area and the average value is obtained by weighting according to the species endangerment level. The phylogenetic diversity conservation index is used to assess the proportion of the total evolutionary history represented by the species set within the target protected area relative to the overall regional population. It is calculated based on a pre-constructed regional species phylogenetic tree and species distribution data in the protected area grid dataset. The region includes the target protected area and the corresponding regions of the non-target protected areas in the same location. The phylogenetic diversity conservation index is the ratio of the sum of the phylogenetic diversity values ​​of each grid in the protected area grid set within the region to the sum of the phylogenetic diversity values ​​of each grid in the protected area grid set and the control grid set within the region.

4. The method for quantitatively evaluating the ecological effectiveness of global protected areas based on time-series remote sensing core variables as described in claim 3, characterized in that, The indicators for the conservation effectiveness dimension include the net change rate of natural vegetation cover, the net change rate of key species populations, and the net reduction rate of human activity disturbance, specifically: The calculation process of the net change rate of natural vegetation cover is as follows: extract the natural vegetation area of ​​the protected area grid set dataset and the control grid set dataset in the base year and the final year, respectively. The natural vegetation area is obtained by statistically analyzing the area of ​​the forest, grassland and wetland types in the natural ecosystem patch data. Calculate the area change rate of each grid. Subtract the area change rate of each grid in the control grid set from the area change rate of each grid in the protected area grid set. The difference is the net change rate. The calculation process for the net population change rate of the key species is as follows: Population monitoring data of the key protected species in the base year and terminal year are obtained from the protected area grid dataset and the control grid dataset. The population monitoring data includes at least one of the relative abundance index monitored by infrared cameras, transect survey encounter rate, or fixed quadrat count. The population change rate of each grid in the protected area grid dataset and each grid in the control grid dataset is calculated from the base year to the terminal year. The difference between the population change rate of each grid in the protected area grid dataset and the population change rate of each grid in the control grid dataset is the net population change rate of the key species. The calculation process of the net reduction rate of human activity interference is as follows: based on the interpretation of the multi-source time-series remote sensing satellite image data, or based on the human activity interference intensity obtained from the human activity interference dataset, the human activity interference intensity of the protected area grid set dataset and the control grid set dataset in the base year and the final year are obtained respectively, the rate of change of interference intensity of each is calculated, and then the rate of change of each grid in the control grid set is subtracted from the rate of change of each grid in the protected area grid set, and the difference is the net reduction rate.

5. The method for quantitatively evaluating the ecological effectiveness of global protected areas based on time-series remote sensing core variables as described in claim 1, characterized in that, The propensity score matching method specifically includes: Based on the geospatial dataset, environmental variables of each grid in the protected area grid set and the control grid set are extracted simultaneously as matching covariates. The environmental variables include at least one of altitude, slope, aspect, distance from road, distance from settlement, annual precipitation and average annual temperature. Using a logistic regression algorithm, with the matching covariate as the independent variable and the group to which the evaluation grid belongs as the dependent variable, the tendency score of each grid in the protected area grid set and the control grid set to belong to the protected area grid set is calculated. The nearest neighbor matching method is used to identify the grid in the control grid set that has the smallest absolute difference in matching tendency scores among all grids in the protected area grid set. The matching results are tested for balance. The standardized mean differences of each matching covariate in the matched protected area grid set and the control grid set are calculated. When the absolute value of all standardized mean differences is less than 0.25, the matching results are deemed to meet the covariate distribution consistency requirements.

6. The method for quantitatively assessing the ecological effectiveness of global protected areas based on time-series remote sensing core variables as described in claim 1, characterized in that, The normalization process uses the Min-Max normalization method to linearly map all original index values ​​to the [0,1] interval. Positive indices use the forward normalization formula, while negative indices use the reverse normalization formula, ensuring that all indices are positive after standardization. The larger the value, the better the protection effect.

7. The method for quantitatively evaluating the ecological effectiveness of global protected areas based on time-series remote sensing core variables as described in claim 1, characterized in that, The construction process of the pre-defined weights, determined by the analytic hierarchy process based on the consistency test of the judgment matrix, specifically includes: The weights are general weights that are pre-constructed based on knowledge in the field of ecological protection and can be reused for different target protected areas; A three-layer hierarchical structure is constructed, comprising a target layer, a criterion layer, and an indicator layer. The target layer is the comprehensive score of the protection effectiveness, the criterion layer includes the protection object dimension and the protection effectiveness dimension, and the indicator layer includes all indicator values ​​of the protection object dimension and the protection effectiveness dimension. Based on the preset 1-9 scale method, a criterion layer judgment matrix and an indicator layer judgment matrix are constructed. The elements of the judgment matrix represent the relative importance values ​​of any two indicators within the same level. Calculate the maximum eigenvalue and corresponding eigenvector of each judgment matrix, and normalize the eigenvectors to obtain the weight vectors of each dimension of the criterion layer and the weight vectors of each index value of the index layer. Calculate the consistency ratio of each judgment matrix. When the consistency ratio is less than 0.1, the weight vector is deemed valid.

8. The method for quantitatively assessing the ecological effectiveness of global protected areas based on time-series remote sensing core variables as described in claim 1, characterized in that, The preset level threshold table includes four levels, specifically: Excellent level, corresponding to a comprehensive score of no less than 85 points for the protection effectiveness; A good rating corresponds to a comprehensive score of protection effectiveness of not less than 70 points and less than 85 points. Medium level, corresponding to a comprehensive score of protection effectiveness of not less than 60 points and less than 70 points; A poor rating corresponds to a comprehensive protection effectiveness score of less than 60 points.

9. The method for quantitatively evaluating the ecological effectiveness of global protected areas based on time-series remote sensing core variables as described in claim 1, characterized in that, The results of the quantitative assessment of ecological effectiveness are presented in at least one of the following visualization formats: Radar charts are used to display the standardized score distribution of each indicator value in the dimensions of protected objects and protection effectiveness of the target protected area, and to identify shortcomings in protection effectiveness. Spatial heat map is used to map the protection effectiveness dimension index values ​​of each grid in the grid set of the target protected area to the spatial range of the target protected area, and to identify the spatial heterogeneity of protection effectiveness; The traffic light warning map is used to mark the changing trends of each indicator value of the target protected area in terms of the protection object dimension and the protection effectiveness dimension from the base period to the final state year in three categories: red light indicates that the indicator value has decreased compared with the base year, yellow light indicates that the change in the indicator value compared with the base year is within a preset threshold, and green light indicates that the indicator value has increased compared with the base year.

10. The method for quantitatively evaluating the ecological effectiveness of global protected areas based on time-series remote sensing core variables as described in claim 1, characterized in that: The multi-source time-series remote sensing satellite imagery data comes from at least one of Sentinel series, Landsat series, or MODIS remote sensing satellite products. The human activity disturbance dataset includes the Global Human Footprint dataset; The geospatial dataset includes at least one of digital elevation models, road network data, settlement distribution data, and climate datasets.

Images