Method for identifying key ecological zones based on ecosystem multifunctionality and stability

Abstract

The application relates to the technical field of ecological remote sensing, in particular to a key ecological area identification method based on the multifunctionality and stability of an ecological system. The method comprises the following steps: obtaining a remote sensing image of a target area, and extracting an ecological function feature set according to the remote sensing image; calculating an ecological multifunctionality index based on the ecological function feature set; extracting an ecological state time sequence by using the remote sensing image; calculating an ecological fluctuation intensity according to the ecological state time sequence; identifying an ecological resilience advantage area according to the ecological fluctuation intensity; evaluating ecological restoration force in the ecological resilience advantage area; calculating a carbon cycle rate based on the ecological restoration force, and identifying a restoration potential partition of the ecological resilience advantage area according to the carbon cycle rate. The application strengthens the comprehensive discrimination ability of the key ecological area based on the ecological remote sensing technology, so that the accuracy of the ecological key area identification and the decision support efficiency are improved.

Description

Technical Field

[0001] This invention relates to the field of ecological remote sensing technology, and in particular to a method for identifying key ecological zones based on the multifunctionality and stability of ecosystems. Background Technology

[0002] In existing ecosystem assessment and ecological zoning studies, high demands are placed on the completeness, spatiotemporal consistency, and comprehensive analytical capabilities of ecosystem state characterization. However, existing technologies still have the following shortcomings: First, existing ecological remote sensing analysis methods mostly focus on the extraction and analysis of single ecological indicators, lacking the ability to integrate and model multiple ecological functional elements such as water distribution, carbon cycle processes, and soil structure characteristics, resulting in a single dimension of ecological function evaluation; Second, existing ecological stability analysis methods are usually based on single-phase or low-frequency time-series remote sensing data for change detection, failing to effectively characterize the continuous fluctuation process of ecological state and disturbance response mechanisms, resulting in insufficient accuracy in the dynamic identification of ecosystem resilience and recovery capacity; Third, existing ecological zoning methods mostly rely on static threshold grading or empirical rules to divide ecological functional areas, making it difficult to achieve accurate identification of key ecological zones and graded assessment of their recovery potential. Summary of the Invention

[0003] Therefore, the present invention needs to provide a method for identifying key ecological zones based on the multifunctionality and stability of ecosystems, in order to solve at least one of the above-mentioned technical problems.

[0004] To achieve the above objectives, a method for identifying key ecological zones based on the multifunctionality and stability of ecosystems includes the following steps:

[0005] Step S1: Acquire remote sensing images of the target area and extract an ecological function feature set based on the remote sensing images; calculate the ecological multifunctionality index based on the ecological function feature set;

[0006] Step S2: Extract ecological state time series from remote sensing images; calculate ecological fluctuation intensity based on ecological state time series; identify areas with ecological resilience advantages based on ecological fluctuation intensity;

[0007] Step S3: Assess ecological resilience in areas with ecological resilience advantages; calculate carbon cycle rate based on ecological resilience, and identify restoration potential zones in areas with ecological resilience advantages based on carbon cycle rate;

[0008] Step S4: Divide the ecological multifunctionality areas according to the ecological multifunctionality index; perform spatial overlay based on the ecological multifunctionality areas and restoration potential zones to identify key ecological zones.

[0009] The beneficial effects of this invention are as follows:

[0010] 1. By constructing a multi-layered processing chain for remote sensing image multispectral band decomposition, ecological function feature set extraction, and ecological multifunctionality index calculation, the transformation of ecosystem representation from single index to multifunctional comprehensive characterization is realized. This avoids the problem that traditional ecological assessment methods rely solely on vegetation index or single cover index, enabling the ecosystem's production function, hydrological regulation function, and carbon cycle function to be comprehensively expressed within a unified framework, thereby improving the completeness and comprehensiveness of ecological multifunctionality assessment.

[0011] 2. By introducing ecological state time series analysis, ecological fluctuation intensity calculation, and ecological resilience advantage area identification mechanism, a quantitative characterization of the long-term dynamic change process of the ecosystem is realized. This avoids the problem that traditional static analysis based on single-phase remote sensing data cannot reflect the ecological disturbance and recovery process. It enables the ecosystem stability to be continuously characterized from the time dimension, and improves the temporal sensitivity and dynamic analysis capability of ecological resilience identification.

[0012] 3. By integrating multiple surface process parameters such as soil aeration porosity, soil moisture buffering capacity, and root growth suitability, and combining ecological resilience and carbon cycle rate to construct a zoning model for restoration potential, a synergistic analysis of ecological process parameters and ecological function evolution is achieved. This avoids the problem of traditional ecological zoning methods relying solely on empirical thresholds or single ecological indicators, enabling the identification results of key ecological zones to simultaneously reflect the intensity of ecological function and the level of restoration capacity, thereby improving the mechanistic and accurate nature of spatial identification. Attached Figure Description

[0013] Other features, objects, and advantages of the invention will become more apparent from the following detailed description of non-limiting embodiments with reference to the accompanying drawings:

[0014] Figure 1 This is a schematic diagram of the steps in the method for identifying key ecological zones based on the multifunctionality and stability of ecosystems in this invention.

[0015] Figure 2 This is a schematic diagram of the carbon fixation flux per unit area in one embodiment;

[0016] Figure 3 This is a schematic diagram of a remote sensing image according to one embodiment;

[0017] The realization of the objective, functional features and advantages of the present invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

[0018] The technical method of this invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of this invention. Based on the embodiments of this invention, all other embodiments obtained by those skilled in the art without inventive effort are within the scope of protection of this invention.

[0019] Furthermore, the accompanying drawings are merely illustrative of the invention and are not necessarily drawn to scale. The same reference numerals in the drawings denote the same or similar parts, and therefore repeated descriptions of them will be omitted. Some block diagrams shown in the drawings are functional entities and do not necessarily correspond to physically or logically independent entities. These functional entities can be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor methods and / or microcontroller methods.

[0020] It should be understood that although the terms "first," "second," etc., may be used herein to describe various units, these units should not be limited by these terms. These terms are used merely to distinguish one unit from another. For example, without departing from the scope of the exemplary embodiments, a first unit may be referred to as a second unit, and similarly, a second unit may be referred to as a first unit. The term "and / or" as used herein includes any and all combinations of one or more of the associated listed items.

[0021] To achieve the above objectives, please refer to Figures 1 to 3 This invention provides a method for identifying key ecological zones based on the multifunctionality and stability of ecosystems, the method comprising the following steps:

[0022] Step S1: Acquire remote sensing images of the target area and extract an ecological function feature set based on the remote sensing images; calculate the ecological multifunctionality index based on the ecological function feature set;

[0023] In one embodiment, multi-temporal remote sensing images (such as Sentinel-2 or Landsat 8) are acquired for the target area, and radiometric calibration and atmospheric correction are performed on the remote sensing images to obtain a corrected remote sensing image dataset. Multispectral band decomposition is performed on the corrected remote sensing images to extract red band, near-infrared band, and shortwave infrared band data. The Normalized Difference Vegetation Index (NDVI) is calculated based on the red band and near-infrared band, and a vegetation cover feature map is generated accordingly. At the same time, a water body index (such as NDWI) is calculated based on the near-infrared band and shortwave infrared band to extract water body distribution features. The vegetation cover features and water body distribution features are spatially fused and standardized to construct an ecological function feature set. The indicators in the ecological function feature set are normalized and weighted to obtain an ecological multifunctionality index.

[0024] In another embodiment, assuming the target area is selected from remote sensing data with an image spatial resolution of 10m, the calculated NDVI range for the red and near-infrared bands is [0.12, 0.86], with an average value of 0.63; the NDWI range for the water body index is... The water coverage rate is approximately 18.5%; the average vegetation coverage is approximately 72%; and the average water distribution index is approximately 0.38. An ecological function vector was constructed after normalization. (corresponding to vegetation function, hydrological function and comprehensive ecological response intensity respectively), the ecological multifunctionality index was calculated to be 0.67 using weighting coefficients of 0.5, 0.3 and 0.2.

[0025] Step S2: Extract ecological state time series from remote sensing images; calculate ecological fluctuation intensity based on ecological state time series; identify areas with ecological resilience advantages based on ecological fluctuation intensity;

[0026] In one embodiment, continuous time-series remote sensing images (e.g., monthly or decadal scale) are extracted from the target area, and an ecological status time series is constructed based on NDVI or a comprehensive ecological index. A sliding window analysis is performed on the time series data to calculate the magnitude of ecological index changes between adjacent time nodes, and the standard deviation, range, and rate of change are statistically analyzed to form an ecological fluctuation intensity parameter. When the ecological fluctuation intensity is below a preset threshold, it is determined to be an ecologically low-disturbance area. Within the ecologically low-disturbance area, analysis is performed in conjunction with surface conditions, including soil pore structure, water retention status, and vegetation stability, thereby identifying areas with ecological resilience advantages.

[0027] In another embodiment, it is assumed that a continuous 12-period NDVI time series is The calculated standard deviation is approximately 0.021, the range is 0.08, and the average rate of change is approximately 0.007 / month. With an ecological fluctuation threshold of 0.03, the region as a whole is in a low-fluctuation state, with approximately 85% of the spatial units meeting the low-disturbance condition. Furthermore, the area with continuous vegetation cover and stability is identified as an area with ecological resilience advantages.

[0028] Step S3: Assess ecological resilience in areas with ecological resilience advantages; calculate carbon cycle rate based on ecological resilience, and identify restoration potential zones in areas with ecological resilience advantages based on carbon cycle rate;

[0029] In one embodiment, vegetation index time series are extracted within ecologically resilient advantage areas to identify the time points of ecological disturbance and continuously track the vegetation recovery process after disturbance; the recovery cycle is calculated based on the time required for the vegetation index to recover to the pre-disturbance level, and the recovery efficiency coefficient is calculated based on the recovery cycle and recovery magnitude to assess ecological resilience; the biomass increment per unit time is estimated based on resilience, and the carbon fixation flux per unit area is calculated in combination with the carbon conversion coefficient, and the net carbon accumulation is calculated in combination with the soil organic carbon background value to obtain the carbon cycle rate; spatial classification is performed based on the carbon cycle rate to delineate areas with recovery potential.

[0030] In another embodiment, it is assumed that three significant disturbance events occur in a region with ecological resilience, with an average recovery period of 18 days and a recovery efficiency coefficient of 0.78; the biomass increment per unit time is approximately If the carbon conversion factor is taken as 0.47, then the carbon fixation flux per unit area is approximately Based on the soil organic carbon baseline value of 2.8%, the net carbon accumulation is approximately [missing information]. The carbon cycle rate was ultimately determined to be 0.42. Based on this, the region with a recovery potential index of approximately 0.81 was classified as a high recovery potential region, the region with an index of 0.55 as a medium recovery potential region, and the region with an index of 0.33 as a low recovery potential region.

[0031] Step S4: Divide the ecological multifunctionality areas according to the ecological multifunctionality index; perform spatial overlay based on the ecological multifunctionality areas and restoration potential zones to identify key ecological zones.

[0032] In one embodiment, the ecological multifunctionality index is standardized to obtain an ecological function spatial distribution map, and the region is divided into spatial units based on the dominant ecological function type to form different functional dominant areas; at the same time, the restoration potential zoning results are rasterized and spatially overlaid with the ecological multifunctional area; the spatial area with high ecological function and high restoration potential is identified by the overlay results as the key ecological area.

[0033] In another embodiment, it is assumed that the average value of the ecological multifunctionality index after standardization is 0.68, of which the production function-dominant area accounts for 42%, the hydrological regulation-dominant area accounts for 35%, and the carbon sink-dominant area accounts for 23%; in the restoration potential zoning, the high potential area accounts for 38%, the medium potential area accounts for 44%, and the low potential area accounts for 18%; after spatial superposition, about 21% of the area simultaneously meets the condition that "ecological function index > 0.7 and restoration potential index > 0.75". This part of the area is identified as the key ecological area, which is mainly distributed in the continuous vegetation area of ​​the river valley and the low-disturbance high-coverage ecological patch area.

[0034] Of particular importance, step S4, which involves delineating ecologically multifunctional zones based on ecological multifunctionality indicators, includes:

[0035] Standardization is performed on the ecological multifunctionality index to obtain a standardized function vector; the contribution of each function index is calculated based on the standardized function vector; and the function index with the largest contribution is extracted based on the contribution of each function index.

[0036] In one embodiment, ecosystem multifunctionality indicators are standardized, including carbon sequestration capacity, soil water retention capacity, vegetation cover, species diversity index, and nutrient cycling efficiency. The Min-Max standardization method is used to map each indicator to the [0,1] interval, resulting in a standardized functional vector. Based on this, a principal component contribution weighting model is introduced to calculate the contribution of each functional indicator: ,in These are the corresponding indicator weighting coefficients. Based on the ranking of each indicator's contribution, the functional indicator with the highest contribution is extracted as the dominant ecological functional indicator.

[0037] In another embodiment, assuming five functional indicators are collected for a certain ecological area, the original data and standardized results are as follows: carbon sequestration capacity: 0.82; soil water retention capacity: 0.65; vegetation cover: 0.74; species diversity index: 0.58; nutrient cycling efficiency: 0.69. The weights are set to 0.25, 0.20, 0.20, 0.15, and 0.20, respectively. The calculated contributions of each indicator are: carbon sequestration capacity: 0.29; soil water retention capacity: 0.15; vegetation cover: 0.20; species diversity index: 0.10; nutrient cycling efficiency: 0.16. The indicator with the largest contribution is carbon sequestration capacity (0.29), which is determined as the dominant ecological function indicator, representing that the dominant ecological characteristic of this area is carbon sequestration.

[0038] The maximum contribution function index is taken as the dominant ecological function type; spatial unit type division is carried out based on the dominant ecological function type to form different functional dominant areas; different functional dominant areas are taken as ecological multifunctional areas.

[0039] In one embodiment, the study area is divided into gridded spatial regions (e.g., 30m×30m or 100m×100m resolution) based on dominant ecological function indicators. Using the dominant function indicator values ​​as the core classification criterion, and combining neighborhood consistency constraints with spatial clustering algorithms (such as K-means or region growth method), spatial units are divided into different functional dominant zones, including carbon sink dominant zones, hydrological regulation dominant zones, biodiversity dominant zones, and nutrient cycling dominant zones.

[0040] In another embodiment, assuming a region is divided into 100 spatial grid units, with carbon sequestration capacity as the maximum contribution indicator, spatial classification is performed according to functional thresholds: 32 grid units with a carbon sequestration capacity > 0.75 form 3 consecutive patches of carbon sink-dominant zones, with a maximum area of ​​approximately 14.2 km²; 41 grid units with a carbon sequestration capacity between 0.5 and 0.75 form 5 patches of transitional ecological functional zones, with an area range of 2.1–9.3 km²; and 27 grid units with a carbon sequestration capacity < 0.5, mainly distributed in construction-affected areas and bare land restoration areas, form 4 patches of low-functional dominant zones, with a maximum area of ​​approximately 3.6 km². These different functional dominant zones collectively constitute a multifunctional ecological region system, where the carbon sink-dominant zone serves as the core ecological functional zone, the transitional zone as a functional buffer zone, and the low-functional zone as a key area for restoration and intervention.

[0041] It should be noted that you should refer to [link / reference]. Figure 2 The horizontal axis represents the east-west direction, and the vertical axis represents the north-south direction. The units of measurement are... The color scale gradually changes from light pink (0.0) to dark green (1.0), with higher values ​​representing stronger carbon fixation flux. Spatially, it exhibits a significant diagonal differentiation feature. The 0-1 interval in the north-south direction, the 0-2 interval in the east-west direction, the 4-6 interval in the north-south direction, and the 3-6 interval in the east-west direction are high-value core areas, dominated by dark green, with carbon fixation flux close to 1.0. The 0-4 interval in the north-south direction, the 3-5 interval in the east-west direction, and some areas in the upper left corner are low-value areas, dominated by light pink, with values ​​close to 0.0. Overall, it shows a significant difference in carbon fixation capacity in different directions.

[0042] It should be noted that you should refer to [link / reference]. Figure 3 The main color is green, with bright green corresponding to vegetation coverage areas, presenting different green gradients, reflecting the differences in surface vegetation density or type; the dark blue area represents water bodies such as rivers and lakes, including a large, irregularly shaped lake in the lower left corner, a winding river in the upper right corner, and a wide body of water in the lower right corner.

[0043] Preferably, in step S1, acquiring remote sensing images of the target area and extracting an ecological function feature set based on the remote sensing images includes:

[0044] Acquire remote sensing images of the target area and decompose the multispectral bands of the remote sensing images; extract the red band, near-infrared band, and short-wave infrared band based on the multispectral bands;

[0045] In one embodiment, multi-temporal remote sensing image data of the target area is acquired, and radiometric calibration and atmospheric correction are performed on the remote sensing images to obtain corrected remote sensing image data. Band separation processing is performed on the corrected remote sensing images to split the original multispectral images into red band (Red), near-infrared band (NIR), and short-wave infrared band (SWIR) data raster. Spatial registration and noise filtering are performed on each band to ensure consistency in spatial resolution and pixel position of different bands, thereby forming a standardized multispectral band data set.

[0046] In another embodiment, assuming the spatial resolution of the acquired remote sensing image of the target area is 10m and the image coverage is approximately 50km², after radiometric calibration, the red band DN value range is [120, 680], the near-infrared band DN value range is [150, 920], and the shortwave infrared band DN value range is [80, 610]. After atmospheric correction, the average red reflectance is approximately 0.32, the average near-infrared reflectance is approximately 0.54, and the average shortwave infrared reflectance is approximately 0.41. After band separation, three sets of raster data are formed, of which approximately [number] effective pixels are... The data consistency registration error is controlled within 0.3 pixels.

[0047] Vegetation indices are calculated based on the red and near-infrared bands; vegetation coverage characteristics are determined using the vegetation indices; water body indices are calculated based on the near-infrared and short-wave infrared bands; and water body distribution characteristics are determined using the water body indices.

[0048] The characteristics of vegetation coverage and water distribution are integrated into a set of ecological function characteristics.

[0049] In one embodiment, the Normalized Difference Vegetation Index (NDVI) is calculated pixel-by-pixel based on the red and near-infrared bands, and the NDVI results are spatially smoothed to reduce noise impact. Vegetation coverage levels are divided according to the NDVI value range to form a vegetation coverage feature map. At the same time, water body indices (such as NDWI or MNDWI) are calculated based on the near-infrared and short-wave infrared bands to identify water bodies and wetland distribution areas. The vegetation coverage features and water body distribution features are spatially superimposed and normalized to construct an ecological function feature set.

[0050] In another embodiment, assuming that the NDVI calculated from red light and near-infrared light ranges from -0.15 to 0.88, with an average value of approximately 0.61, high vegetation cover areas (NDVI > 0.7) account for approximately 46%, medium cover areas account for approximately 38%, and low cover areas account for approximately 16%; the water body index NDWI ranges from -0.52 to 0.63, with water body pixels accounting for approximately 19.3% and wetland transition areas accounting for approximately 11.5%; after normalization, the mean value of vegetation cover features is approximately 0.68, and the mean value of water body distribution features is approximately 0.42; finally, these are fused to form an ecological function feature vector. .

[0051] Preferably, acquiring remote sensing images of the target area and decomposing the multispectral bands of the remote sensing images includes:

[0052] Acquire remote sensing images of the target area, and perform radiometric calibration and atmospheric correction on the remote sensing images to obtain corrected remote sensing images; perform band separation on the corrected remote sensing images to obtain band data;

[0053] In one embodiment, multi-temporal remote sensing image data of the target area is acquired, and radiometric calibration processing is performed on the original remote sensing images to convert DN values ​​into surface reflectance data. An atmospheric correction model (such as the 6S or FLAASH model) is used to eliminate the effects of aerosol scattering and water vapor absorption, resulting in a highly consistent corrected remote sensing image. The corrected remote sensing image is then subjected to multispectral band separation processing, splitting the image data into independent band data such as red band, green band, near-infrared band, and shortwave infrared band. Radiometric consistency normalization processing is then performed on each band to form a standardized band dataset.

[0054] In another embodiment, assuming the acquired remote sensing image has a spatial resolution of 10m and a coverage area of ​​approximately 60km², the original red band DN value range is [110, 720], and the near-infrared band DN value range is [130, 950]. After radiometric calibration, the average red reflectance is approximately 0.28, and the average near-infrared reflectance is approximately 0.56. After atmospheric correction, the aerosol effect is reduced by approximately 35%, and the water vapor absorption error is controlled within 0.02. After band separation, three types of band data gratings are obtained: red, near-infrared, and short-wave infrared, with approximately [number of effective pixels]. The inter-band registration error is less than 0.25 pixels.

[0055] Perform region clipping and spatial resampling based on band data to generate multispectral bands.

[0056] In one embodiment, target region boundary vector information is introduced based on band data, and spatial clipping is performed on the multispectral band data to extract an effective subset of bands within the target region. Based on the original spatial resolution and the target analysis scale requirements, a resampling ratio coefficient is calculated, and bilinear interpolation or nearest neighbor interpolation methods are used to spatially resample the band subsets to unify the spatial resolution of different bands. Raster alignment is then performed on the resampled band data to ensure strict consistency of each band in the spatial coordinate system, generating a standard multispectral band dataset that can be used for multi-source fusion analysis.

[0057] In another embodiment, assuming the original remote sensing image has a spatial resolution of 10m, the target analysis resolution is set to 5m, and the resampling ratio is 2; the effective area after cropping the target region boundary is approximately 48km², corresponding to approximately 1.92 × 10⁻⁶ raster pixels. 6 The spatial alignment error of each band was controlled within 0.3 pixels after the bilinear interpolation method was used during the resampling process. After resampling, the mean value of the red band was adjusted from 0.31 to 0.30, the mean value of the near-infrared band was adjusted from 0.58 to 0.57, and the mean value of the short-wave infrared band was adjusted from 0.42 to 0.41.

[0058] Preferably, performing region clipping and spatial resampling based on band data to generate multispectral bands includes:

[0059] Obtain the target region boundary vector data; perform spatial clipping on the band data based on the target region boundary vector data to obtain a subset of the target region bands; determine the original spatial resolution based on the subset of the target region bands;

[0060] In one embodiment, vector boundary data of the target area (such as Shapefile or GeoJSON format) is acquired and spatially mapped to the multispectral band data using a coordinate system (such as WGS84 or UTM projection). Based on the boundary vector, spatial clipping is performed on the original band raster data, retaining only the valid image data within the boundary range to form a subset of bands in the target area. The spatial resolution parameters of the clipped band subset are calculated, and the original spatial resolution is determined by resolving the geographic spacing of raster pixels (such as pixel size and row / column spacing). The resolution consistency between different bands is recorded to provide basic parameters for subsequent scale unification processing.

[0061] In another embodiment, it is assumed that the target area boundary covers an area of ​​approximately 52 km², corresponding to approximately 350 vector boundary nodes; the original remote sensing image has a spatial resolution of 10 m, and the number of effective pixels after cropping is approximately... The original spatial resolution was calculated to be 10.02m through pixel geographic spacing (with an interpolation error of 0.2%). The maximum spatial offset error for different bands was about 0.3 pixels, mainly concentrated in areas with complex boundaries (such as the junction of rivers and forests), but the overall cropping results maintained spatial consistency.

[0062] The resampling ratio coefficient is calculated based on the original spatial resolution and the preset resolution parameters; spatial resampling is performed on a subset of bands in the target area using the resampling ratio coefficient to obtain resampled bands; grid alignment is performed based on the resampled bands to construct multispectral bands.

[0063] In one embodiment, a scale transformation relationship is established based on the original spatial resolution and the preset target resolution parameters, a resampling ratio coefficient is calculated, and spatial resampling processing is performed on a subset of bands in the target region based on the coefficient. During the resampling process, bilinear interpolation or cubic convolution interpolation methods are used to perform a unified scale transformation on different bands to make the grid size of each band consistent. Grid alignment processing is performed on the resampled multi-band data, including spatial reference unification, pixel grid alignment, and overlapping area correction, thereby constructing a spatially consistent multispectral band dataset.

[0064] In another embodiment, assuming the original spatial resolution is 10m and the preset target resolution is 5m, the resampling ratio is 2; the number of pixels in the cropped band subset is approximately The number of pixels increased to approximately [number] after resampling. During the resampling process, bilinear interpolation was used to control the spatial error of each band within 0.25 pixels. After grid alignment, the maximum spatial offset between the red band and the near-infrared band was reduced from 0.4 pixels to 0.1 pixels, and the overall spatial consistency was improved by about 75%, ensuring that the subsequent calculation of vegetation index and water index has the same spatial reference.

[0065] Preferably, step S2, which identifies areas with ecological resilience advantages based on the intensity of ecological fluctuations, includes:

[0066] If the intensity of ecological fluctuations is less than the preset ecological fluctuation threshold, it is determined to be an ecologically low disturbance area; soil aeration porosity is detected in the ecologically low disturbance area; root growth suitability is assessed based on soil aeration porosity.

[0067] In one embodiment, an ecological fluctuation intensity index is calculated based on the ecological state time series of remote sensing images and compared with a preset ecological fluctuation threshold. When the ecological fluctuation intensity is lower than the threshold, the target spatial unit is determined to be an ecologically low-disturbance area. Within this low-disturbance area, a soil structure parameter inversion model is introduced to detect soil aeration porosity. Specifically, the total soil porosity is calculated by analyzing the proportion of soil particles (sand, silt, and clay) and their spatial distribution characteristics, and the effective aeration porosity is obtained by separating the non-capillary pore portion. Based on this, a root growth suitability evaluation function is constructed, taking aeration porosity, water and air exchange capacity, and pore connectivity as input variables to comprehensively evaluate the plant root penetration capacity and respiratory environment suitability, thereby outputting a root growth suitability index.

[0068] In another embodiment, assuming the ecological fluctuation intensity of a certain ecological unit is 0.32 and the preset threshold is 0.5, it is determined to be an ecologically low disturbance area. Within this area, the content of sand particles is detected to be 38%, silt particles 42%, and clay particles 20%, corresponding to a loam soil type. The total soil porosity is calculated to be approximately 52%, of which the effective aeration porosity is approximately 18%–22%. Combining the pore connectivity index of 0.76 and the oxygen diffusion coefficient of 0.68, the root growth suitability is calculated to be 0.81 (out of 1.0), indicating that the area has relatively good root expansion and respiration environment conditions.

[0069] Soil moisture content is measured in ecologically low-disturbance areas; soil moisture buffering capacity is assessed based on soil moisture content; and ecologically resilient areas are identified based on soil moisture buffering capacity and root growth suitability.

[0070] In one embodiment, soil moisture time series is calculated jointly based on multi-temporal remote sensing inversion data and ground soil moisture model in ecologically low disturbance areas. By statistically analyzing the fluctuation range of moisture content, the frequency of wet-dry cycles, and the moisture retention capacity, a dynamic stability index of moisture is constructed. Moisture fluctuation characteristics are calculated by the standard deviation, range, and rate of change of moisture content, and a moisture stability rate model is introduced to evaluate the soil moisture buffering capacity, that is, the ability of soil to maintain moisture balance and restore the original moisture state under external climate disturbances. The moisture buffering capacity is coupled with root growth suitability for analysis, and a weighted fusion model is used to identify areas with ecological resilience advantages, so as to achieve a comprehensive judgment on ecological stability and vegetation restoration capacity.

[0071] In another embodiment, assuming the soil moisture content of a certain ecologically low-disturbance area at 10 time points is [missing information]. The calculated standard deviation of moisture content was 1.35%, the range was 4%, the frequency of dry-wet transition was 2 times / cycle, and the average rate of moisture change was approximately 0.42% / day; thus, the moisture stability rate was 0.78; combined with the aforementioned root growth suitability of 0.81, the ecological resilience index was calculated to be 0.795 using the comprehensive evaluation function (0.5 × suitability + 0.5 × moisture stability rate); when the ecological resilience index is greater than 0.75, the area is identified as an area with ecological resilience advantage, indicating that the area has both good water regulation capacity and root growth support capacity.

[0072] Preferably, the detection of soil aeration porosity in ecologically low-disturbance areas includes:

[0073] The proportion of soil particle composition was measured in ecologically low-disturbance areas, including the content of sand, silt, and clay; the soil texture type was determined based on the proportion of soil particle composition.

[0074] In one embodiment, within a determined low-disturbance ecological area, soil particle composition parameters, including the mass or volume fractions of sand, silt, and clay, are obtained by collecting soil samples or by inverting soil spectral characteristics based on remote sensing. The proportions of the three particle types are used as input features, and the soil texture type is determined based on the soil triangulation method (USDA soil texture classification system), thereby classifying the soil into types such as sandy soil, loam, or clay. The texture type within the area is then rasterized using spatial interpolation methods to form a continuously distributed spatial layer of soil texture.

[0075] In another embodiment, assuming that the soil particle composition ratio in a certain ecologically low-disturbance area is: 35% sand, 45% silt, and 20% clay; the soil texture type of the area is determined to be "loam" based on the soil triangulation diagram; the spatial scale statistics show that: loam accounts for 72%, sandy loam accounts for 18%, and clay loam accounts for 10%; among which, the loam area is mainly distributed in gentle slopes with a slope of 2°–5°, the particle distribution is relatively uniform, and the structural stability is relatively high.

[0076] The total soil porosity is calculated based on the soil texture type; the total soil porosity is corrected based on the non-capillary pore volume fraction corresponding to the soil texture type, and recorded as soil aeration porosity.

[0077] In one embodiment, the total porosity of the soil is calculated based on empirical parameter models corresponding to different soil texture types. The calculation is based on the soil particle packing density and pore formation law. The total porosity is corrected according to the non-capillary pore volume fraction (i.e., the proportion of aeration pores) corresponding to the soil texture type to eliminate the influence of capillary pores on aeration capacity, thereby obtaining the soil aeration porosity that characterizes the soil's gas exchange capacity. This parameter is used to characterize oxygen diffusion capacity, root respiration conditions, and microbial activity space, providing a key physical indicator for ecological resilience assessment.

[0078] In another embodiment, assuming that the empirical model for loam soil calculates the total soil porosity to be 48%, with the non-capillary pore volume fraction ranging from 18% to 22%, and taking the median of 20% as a correction coefficient, the corrected soil aeration porosity is 48% × 0.20 = 9.6%. Spatially, the aeration porosity range is 8.5%–11.2%, with approximately 8.7% in the top of the slope, approximately 10.3% in the middle of the slope, and approximately 11.1% in the foot of the slope. The results indicate that the aeration conditions are better in the foot of the slope, which is conducive to root expansion and oxygen exchange.

[0079] Preferably, assessing soil moisture buffering capacity based on soil moisture content includes:

[0080] The standard deviation and range of soil moisture content are calculated based on the soil moisture content, and these are used as parameters for moisture fluctuation amplitude. The rate of change of moisture content between adjacent time points is calculated based on the parameters for moisture fluctuation amplitude.

[0081] In one embodiment, firstly, based on multi-temporal soil moisture monitoring data in an ecologically low-disturbance area, a moisture content time series is constructed. Statistical analysis is performed on this time series to calculate the standard deviation and range of moisture content, which characterize the dispersion and fluctuation range of soil moisture in the time dimension. The standard deviation and range are then used as parameters for the amplitude of moisture fluctuation. Based on the difference in moisture content and the time interval between adjacent time nodes, the rate of change of moisture content is calculated to characterize the intensity of dynamic changes in soil moisture on a short time scale, thereby forming a quantitative description system for the process of moisture fluctuation and change.

[0082] In another embodiment, assuming that the soil moisture content of a certain ecologically low-disturbance area is [18.2%, 19.0%, 18.5%, 20.1%, 17.8%, 18.3%, 19.5%, 18.9%, 18.1%, 20.0%] for 10 consecutive time points, the standard deviation of the moisture content is calculated to be 0.82%, and the range is 2.3%. After using these two as parameters for the amplitude of moisture fluctuation, the difference calculation is performed on adjacent time points to obtain the rate of change of moisture content as [0.8%, -0.5%, 1.6%, -2.3%, 0.5%, 1.2%, -0.6%, -0.8%, 1.9%] / time unit, where the maximum rate of change is 2.3% and the average rate of change is 0.99% / time unit, indicating that the moisture fluctuation in this area is relatively gentle overall but there are local short-term disturbances.

[0083] The frequency of wet-dry transitions is counted based on the rate of change in moisture content; the moisture stability rate is calculated based on the frequency of wet-dry transitions and the magnitude of moisture fluctuations; and the soil moisture buffering capacity is assessed based on the moisture stability rate.

[0084] In one embodiment, after obtaining the sequence of water content change rates, a dry-wet state determination rule is set. When the water content changes from above the field capacity threshold to below the threshold, it is determined as a dry-wet transition event, and the number of transitions in the entire time series is counted. Combining the water fluctuation amplitude parameters (standard deviation and range), a water stability calculation model is constructed. High fluctuation corresponds to low stability, and low fluctuation corresponds to high stability, and a normalized mapping is performed to obtain the water stability index. Finally, the water stability index characterizes the strength of the soil's ability to maintain water balance under external climate disturbances, and the soil water buffering capacity, that is, the soil's ability to regulate and buffer precipitation input and evapotranspiration loss, is evaluated accordingly.

[0085] In another embodiment, it is assumed that the moisture content crosses the field capacity threshold (set to 19%) twice during the aforementioned 10 time points, resulting in a wet-dry transition. Combining the moisture fluctuation amplitude parameters (standard deviation 0.82%, range 2.3%), the moisture stability rate is calculated to be 0.79 (range 0–1). Spatially, the moisture stability rates of different sub-regions are 0.73 at the top of the slope, 0.81 in the middle of the slope, and 0.85 at the foot of the slope. According to the set rules, when the moisture stability rate is greater than 0.75, it is determined to be a region with strong moisture buffering capacity. Therefore, the middle and foot of the slope regions show better moisture buffering capacity, while the top of the slope region is relatively weak.

[0086] Preferably, step S3, which assesses ecological resilience in areas with ecological resilience advantages, includes:

[0087] Obtain the time series of vegetation index in areas with ecological resilience advantages; identify the time points of index decline based on the vegetation index time series, and use the time points of index decline as the time points of ecological disturbance;

[0088] In one embodiment, the Normalized Difference Vegetation Index (NDVI) is calculated for the target ecological resilience advantage area based on multi-temporal remote sensing images (such as Sentinel-2 or high-resolution UAV images) to form a time series. .right Perform sliding window difference analysis (window length of 5 time steps) to calculate .when When the value falls below a preset threshold (e.g., -0.08), the time point is marked as an ecological disturbance time point. The disturbance points are then clustered and merged (time intervals of less than 3 days are considered as the same disturbance event) to obtain an ecological disturbance event sequence.

[0089] In another embodiment, assuming that a continuous 12-month NDVI time series of a region with ecological resilience advantages is obtained,

[0090] .calculate The subsequent decrease point was: in the 5th month (0.80→0.74). ), 6th month (0.74→0.72, ), where the comprehensive disturbance judgment threshold is If the fifth month is considered the time point for ecological disturbance, then this time point corresponds to a significant disturbance event, such as a short-term drought or localized anthropogenic disturbance.

[0091] Record the changes in vegetation index after the ecological disturbance time point; calculate the vegetation recovery time cycle based on the changes in vegetation index; calculate the recovery efficiency coefficient based on the vegetation recovery time cycle; and assess the ecological resilience based on the recovery efficiency coefficient.

[0092] In one embodiment, starting from the ecological disturbance time point t0, the subsequent vegetation index sequence is extracted. Calculate the recovery increment .when Record the recovery time period when the value recovers to 95% of the baseline value before the disturbance. The recovery efficiency coefficient is defined as follows: ,in To restore the maximum vegetation index within the window. According to Classify ecological resilience: High resilience is defined as 0.015–0.03, medium resilience as <0.015, and low resilience as <0.015.

[0093] In another embodiment, it is assumed that the ecological disturbance occurs in the 5th month. Pre-disturbance reference The value was 0.80. The subsequent recovery sequence was as follows: 0.72 in month 6, 0.73 in month 7, 0.78 in month 8, 0.79 in month 9, and 0.81 in month 10. Therefore, when... The first recovery to ≥0.80 occurred in the 10th month, therefore the recovery time period... Months. The longest recovery period. Calculate the recovery efficiency coefficient: Based on this, it is determined that the restoration efficiency of this area is low, indicating that the ecosystem responds slowly to disturbances, which may be related to soil moisture limitation or monoculture of vegetation. The overall ecological resilience is at a low to medium level.

[0094] Preferably, step S3, which calculates the carbon cycle rate based on ecological resilience and identifies the restoration potential zones of ecologically resilient areas based on the carbon cycle rate, includes:

[0095] Calculate biomass increment per unit time based on ecological resilience; calculate carbon fixation flux per unit area based on biomass increment per unit time; obtain soil organic carbon baseline value;

[0096] In one embodiment, based on the ecological resilience index (Based on the aforementioned vegetation restoration efficiency or) (Obtained from the recovery model), and the biomass increment per unit time is established through empirical models or biomass remote sensing inversion models. For example, using ,in Biomass conversion factor (e.g.) Combined with the photosynthetic carbon conversion efficiency coefficient (e.g., 0.45), calculate the carbon fixation flux per unit area. This allows us to obtain the carbon uptake intensity. Soil organic carbon baseline values ​​are obtained through historical soil monitoring data or global soil databases. .

[0097] In another embodiment, assuming a region with ecological resilience advantage has a resilience index of [missing information] within one year. If it is 0.6, then: Carbon fixation flux Simultaneously, based on regional soil sampling data, The distribution is as follows (unit: ): [85,88,90,87,92], with an average value of approximately This benchmark is used to measure the background level of ecosystem carbon storage and provides an initial reference for carbon gain assessment.

[0098] Of particular importance is that the calculation of carbon fixation flux per unit area based on biomass increment per unit time includes:

[0099] The biomass carbon conversion coefficient is determined based on the biomass increment per unit time; the biomass carbon conversion coefficient and the biomass increment per unit time are multiplied to obtain the carbon fixation per unit time.

[0100] In one embodiment, the dynamic changes in biomass during the vegetation growth process of the ecosystem are continuously monitored to obtain the biomass increment per unit time. (For example, with) Characterization), and combined with the pre-defined biomass carbon conversion coefficient in the vegetation type database. (This represents the proportionality coefficient of conversion from unit biomass to carbon storage). The biomass-carbon conversion coefficient is adjusted based on vegetation type, temperature conditions, and photosynthetic efficiency to obtain the dynamic conversion coefficient. Through the formula: ; Calculate carbon fixation per unit time .

[0101] In another embodiment, assume that the biomass increment per unit time in a certain area is as follows ( The corresponding biomass carbon conversion coefficients for grids A (18.5), B (22.0), C (15.3), and D (27.8) are as follows: ; ; ; The calculated result of carbon fixation per unit time is: ; ; ; This yields the distribution of carbon fixation capacity per unit time in different spatial units.

[0102] Obtain the area of ​​the target spatial unit; normalize the carbon fixation per unit time based on the area of ​​the target spatial unit to obtain the carbon fixation flux per unit area.

[0103] In one embodiment, the area A (unit: m²) of the target spatial unit is obtained, and the carbon fixation per unit time is calculated. The area was normalized to obtain the carbon fixation flux per unit area. : This indicator is used to eliminate spatial scale differences, making the carbon sequestration capacity of different spatial units comparable, and thus used for ecological resilience assessment and spatial zoning analysis.

[0104] In another embodiment, it is assumed that the target area is divided into 4 spatial units with areas of: A1 = 900 m²; A2 = 1200 m²; A3 = 800 m²; A4 = 1500 m². Combined with the carbon fixation per unit time: ; ; ; The carbon fixation flux per unit area are as follows: ; ; ; This yields the distribution of carbon fixation flux per unit area in different spatial units, which can be used to identify regions with high carbon sink potential and regions with low carbon sink sensitivity.

[0105] Net carbon accumulation is calculated based on soil organic carbon baseline and carbon fixation flux per unit area; carbon cycle rate is calculated based on net carbon accumulation; and restoration potential zones are identified based on carbon cycle rate for areas with ecological resilience advantages.

[0106] In one embodiment, net carbon accumulation is defined as ,in The time scale is (year). The soil carbon release coefficient is given (e.g., 0.02). Based on this, the carbon cycle rate is defined as... This is used to characterize the carbon input-output balance efficiency. According to... Size is used to categorize ecological restoration potential, for example: The range is 0.2–0.5 for high recovery potential, 0.2–0.5 for medium recovery potential, and <0.2 for low recovery potential.

[0107] In another embodiment, it is assumed that Year, , , Then: carbon loss term Net carbon accumulation Carbon cycle rate Spatially, assuming three sub-regions, the calculation results are as follows: Region A: (High Recovery Potential Zone) Area B: (Mid-recovery potential zone) area : (Low recovery potential area); this can be used to identify areas. It has strong carbon sequestration capacity and ecological resilience restoration advantages, regional It is in a state of dynamic equilibrium, while the region There is a trend of carbon deficit and weak recovery potential, requiring key ecological restoration interventions.

[0108] Preferably, the zoning for identifying ecologically resilient areas based on carbon cycle rates includes:

[0109] A recovery potential index greater than 0.75 is identified as a high recovery potential region; a recovery potential index between 0.4 and 0.75 is identified as a medium recovery potential region; and a recovery potential index less than 0.4 is identified as a low recovery potential region.

[0110] In one embodiment, the ecosystem resilience potential index is defined as a comprehensive indicator based on the vegetation index change rate, soil moisture response coefficient, and normalized carbon sequestration flux. The ecological resilience advantage area is then divided into grids (e.g., 30m × 30m resolution). The resilience potential index is calculated for each grid cell. The calculation method is as follows: Recovery slope + 0.3 × normalized carbon fixation flux + 0.3 × moisture response stability coefficient.

[0111] Based on the calculation results, the study area is classified and determined: when When it is determined to be a high recovery potential area; when When it is determined to be an area with medium recovery potential; when Areas identified as having low recovery potential are thus identified, thereby creating a spatial distribution map of recovery potential.

[0112] In another embodiment, assuming a certain ecological resilience advantage area with a total of 100 grid cells is selected, the following distribution is obtained after calculating its recovery potential index: There are a total of 28 grid cells; There are a total of 46 grid cells; There are a total of 26 grid cells. Then... The 28 grids were classified as high restoration potential areas, mainly concentrated in riparian buffer zones and areas with continuous vegetation cover; the 46 grids in the 0.4–0.75 range were classified as medium restoration potential areas, mainly distributed in agroforestry ecotones and slightly disturbed areas; The 26 grids are divided into low restoration potential areas, mainly located in bare land and degraded grassland patches.

[0113] The regions with high recovery potential, medium recovery potential, and low recovery potential are designated as recovery potential zones.

[0114] In one embodiment, the three types of regions are spatially clustered and their boundaries are merged to form continuous recovery potential partition units. A region merging algorithm based on adjacency is used to merge adjacent and similar grid units to form three types of spatial patches: high recovery potential partitions, medium recovery potential partitions, and low recovery potential partitions. The area ratio, boundary complexity, and spatial connectivity index of each partition are calculated.

[0115] In another embodiment, it is assumed that the following restoration potential zoning results are obtained after spatial clustering: the high restoration potential zone consists of 6 continuous patches, with the largest patch area being approximately 12.5 km² and the average patch area being approximately 4.3 km²; the medium restoration potential zone consists of 9 patches, with areas ranging from 1.2 km² to 8.7 km²; and the low restoration potential zone consists of 5 patches, mainly fragmented distribution areas, with the largest patch being only 3.1 km². The high restoration potential zone accounts for approximately 32% of the total area, the medium restoration potential zone accounts for 48%, and the low restoration potential zone accounts for 20%. These three types of restoration potential zones together constitute the spatial structure model of ecological resilience, where the high restoration potential zone can serve as the core area for priority ecological restoration and natural restoration, the medium restoration potential zone as the transitional area for auxiliary restoration, and the low restoration potential zone as the key area for governance and intervention.

[0116] Therefore, the embodiments should be considered as exemplary and non-limiting in all respects, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of the equivalents of the application are intended to be included within the invention.

[0117] The above description is merely a specific embodiment of the present invention, enabling those skilled in the art to understand or implement the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the present invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features of the invention herein.

Claims

1. A method for identifying key ecological zones based on the multifunctionality and stability of ecosystems, characterized in that, Includes the following steps: Step S1: Acquire remote sensing images of the target area and extract an ecological function feature set based on the remote sensing images; calculate the ecological multifunctionality index based on the ecological function feature set; Step S2: Extract ecological state time series from remote sensing images; calculate ecological fluctuation intensity based on ecological state time series; identify areas with ecological resilience advantages based on ecological fluctuation intensity; Step S3: Assess ecological resilience in areas with ecological resilience advantages; calculate carbon cycle rate based on ecological resilience, and identify restoration potential zones in areas with ecological resilience advantages based on carbon cycle rate; Step S4: Divide the ecological multifunctional areas according to the ecological multifunctionality index; Spatial overlay is performed based on ecological multifunctional zones and restoration potential zones to identify key ecological areas.

2. The method for identifying key ecological zones based on the multifunctionality and stability of ecosystems according to claim 1, characterized in that, Step S1 involves acquiring remote sensing images of the target area and extracting an ecological function feature set based on the remote sensing images, including: Acquire remote sensing images of the target area and decompose the multispectral bands of the remote sensing images; extract the red band, near-infrared band, and short-wave infrared band based on the multispectral bands; Vegetation indices are calculated based on the red and near-infrared bands; vegetation coverage characteristics are determined using the vegetation indices; water body indices are calculated based on the near-infrared and short-wave infrared bands; and water body distribution characteristics are determined using the water body indices. The characteristics of vegetation coverage and water distribution are integrated into a set of ecological function characteristics.

3. The method for identifying key ecological zones based on the multifunctionality and stability of ecosystems according to claim 2, characterized in that, Acquiring remote sensing images of the target area and decomposing the multispectral bands of the remote sensing images includes: Acquire remote sensing images of the target area, and perform radiometric calibration and atmospheric correction on the remote sensing images to obtain corrected remote sensing images; perform band separation on the corrected remote sensing images to obtain band data; Perform region clipping and spatial resampling based on band data to generate multispectral bands.

4. The method for identifying key ecological zones based on the multifunctionality and stability of ecosystems according to claim 3, characterized in that, Performing region clipping and spatial resampling based on band data to generate multispectral bands includes: Obtain the target region boundary vector data; perform spatial clipping on the band data based on the target region boundary vector data to obtain a subset of the target region bands; determine the original spatial resolution based on the subset of the target region bands; The resampling ratio coefficient is calculated based on the original spatial resolution and the preset resolution parameters; spatial resampling is performed on a subset of bands in the target area using the resampling ratio coefficient to obtain resampled bands; grid alignment is performed based on the resampled bands to construct multispectral bands.

5. The method for identifying key ecological zones based on the multifunctionality and stability of ecosystems according to claim 1, characterized in that, Step S2, which identifies areas with ecological resilience advantages based on the intensity of ecological fluctuations, includes: If the intensity of ecological fluctuations is less than the preset ecological fluctuation threshold, it is determined to be an ecologically low disturbance area; soil aeration porosity is detected in the ecologically low disturbance area; root growth suitability is assessed based on soil aeration porosity. Soil moisture content is measured in ecologically low-disturbance areas; soil moisture buffering capacity is assessed based on soil moisture content; and ecologically resilient areas are identified based on soil moisture buffering capacity and root growth suitability.

6. The method for identifying key ecological zones based on the multifunctionality and stability of ecosystems according to claim 5, characterized in that, Testing soil aeration porosity in ecologically low-disturbance areas includes: The proportion of soil particle composition was measured in ecologically low-disturbance areas, including the content of sand, silt, and clay; the soil texture type was determined based on the proportion of soil particle composition. The total soil porosity is calculated based on the soil texture type; the total soil porosity is corrected based on the non-capillary pore volume fraction corresponding to the soil texture type, and recorded as soil aeration porosity.

7. The method for identifying key ecological zones based on the multifunctionality and stability of ecosystems according to claim 5, characterized in that, Soil moisture buffering capacity is assessed based on soil moisture content, including: The standard deviation and range of soil moisture content are calculated based on the soil moisture content, and these are used as parameters for moisture fluctuation amplitude. The rate of change of moisture content between adjacent time points is calculated based on the parameters for moisture fluctuation amplitude. The frequency of wet-dry transitions is counted based on the rate of change in moisture content; the moisture stability rate is calculated based on the frequency of wet-dry transitions and the magnitude of moisture fluctuations; and the soil moisture buffering capacity is assessed based on the moisture stability rate.

8. The method for identifying key ecological zones based on the multifunctionality and stability of ecosystems according to claim 1, characterized in that, Step S3, which assesses ecological resilience in areas with ecological resilience advantages, includes: Obtain the time series of vegetation index in areas with ecological resilience advantages; identify the time points of index decline based on the vegetation index time series, and use the time points of index decline as the time points of ecological disturbance; Record the changes in vegetation index after the ecological disturbance time point; calculate the vegetation recovery time cycle based on the changes in vegetation index; calculate the recovery efficiency coefficient based on the vegetation recovery time cycle; and assess the ecological resilience based on the recovery efficiency coefficient.

9. The method for identifying key ecological zones based on the multifunctionality and stability of ecosystems according to claim 1, characterized in that, Step S3 involves calculating the carbon cycle rate based on ecological resilience and identifying the restoration potential zones of ecologically resilient areas based on the carbon cycle rate, including: Calculate biomass increment per unit time based on ecological resilience; calculate carbon fixation flux per unit area based on biomass increment per unit time; obtain soil organic carbon baseline value; Net carbon accumulation is calculated based on soil organic carbon baseline and carbon fixation flux per unit area; carbon cycle rate is calculated based on net carbon accumulation; and restoration potential zones are identified based on carbon cycle rate for areas with ecological resilience advantages.

10. The method for identifying key ecological zones based on the multifunctionality and stability of ecosystems according to claim 9, characterized in that, Based on carbon cycle rates, regions with ecological resilience advantages and restoration potential are identified as follows: A recovery potential index greater than 0.75 indicates a high recovery potential region; a recovery potential index between 0.4 and 0.75 indicates a medium recovery potential region; and a recovery potential index below 0.4 indicates a low recovery potential region. The regions with high recovery potential, medium recovery potential, and low recovery potential are designated as recovery potential zones.

Images