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

Abstract

The application discloses a key ecological zone identification method based on ecosystem multifunctionality and stability. The method combines the multifunctionality and stability of an ecosystem to determine the key ecological zone pattern of a region, and inputs limited resources and funds into the key ecological zone of the region, so that the maximization of the ecological protection benefit of the region is realized.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of regional ecosystem management and regional land use planning and management technology, and particularly relates to a key ecological area identification method based on the multi-functionality and stability of an ecological system. BACKGROUND

[0002] Key ecological areas (CEAs) are important areas that maintain biodiversity and ecosystem functions, and continuously supply ecosystem services to human society, playing an important role in regional ecological security. However, in recent years, under the influence of climate change and human activities, the function of regional key ecological areas has been deteriorating, and the problem of declining supply capacity of ecosystem services cannot be ignored. In the process of regional ecological protection policy making and restoration, due to the limitation of investment funds and resources, determining the protection order, identifying regional key ecological areas and priority protection hotspots are the most basic and critical stages. Investing limited resources and funds into regional key ecological areas is conducive to improving the efficiency of ecological protection and management, so the identification of regional key ecological areas has become an important tool for improving ecosystem degradation and maintaining reasonable landscape planning, and the related research has been widely concerned in the global range. Therefore, in-depth research on the identification of regional key ecological areas will provide a fundamental solution for coordinating regional ecological protection and sustainable social and economic development.

[0003] Ecosystem services are the bridge connecting the natural system and human well-being. Ecosystems continuously provide products and services for human beings through ecological functions to meet the needs of human well-being. At present, the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services (IPBES) in the latest survey report pointed out that the ability of nature to maintain good quality of life has declined significantly worldwide since 1970. Therefore, it is of great practical significance to carry out research on the identification of key ecological areas based on ecosystem services. Key ecological areas based on ecosystem services are defined as "areas with important ecological service function supply compared to other areas", which should be fully protected. From the perspective of ecosystem services, the identification of regional key ecological areas is increasingly valued, and the research results will provide important technical support for the improvement and efficiency of regional ecosystem service supply.

[0004] At present, in the field of key ecological zone identification, most of them are based on the static evaluation of elements, such as 3S technology and model simulation methods, and the spatial and temporal superposition analysis of various single ecological service functions, focusing on the static evaluation or long-term average evaluation of ecological system service supply. Such methods pay too much attention to the static or average characteristics of the ecological system, and do not consider the dynamic change characteristics of the long-term ecological system. As an important intrinsic characteristic of the long-term ecological system, the dynamic stability and multifunctionality of the long-term ecological system are beneficial to realize the stable and efficient supply of regional ecological services, but the related research has not systematically introduced the concepts and determination methods of ecological system multifunctionality and stability into the framework of regional key ecological zone identification, which needs to be further deepened.

[0005] Based on the above problems, the present inventors have made in-depth research on the dynamic stability and multifunctionality of the long-term ecological system, in the hope of discovering a method for quickly and accurately obtaining the key ecological zone. SUMMARY

[0006] In order to overcome the above problems, the present inventors have made intensive research and designed a key ecological zone identification method based on the multifunctionality and stability of the ecological system. In this method, the key ecological zone pattern of the region is determined from the two aspects of the multifunctionality and stability of the ecological system, and the limited resources and funds are invested in the key ecological zone of the region, so as to maximize the regional ecological protection benefit, thereby completing the present application.

[0007] Specifically, the purpose of the present application is to provide a key ecological zone identification method based on the multifunctionality and stability of the ecological system, which comprises the following steps:

[0008] Step 1: Obtain a multi-source remote sensing data set of the target region within N years;

[0009] Step 2: Quantize the spatial distribution pattern of each ecological system service in the target region within N years, and extract the spatial distribution pattern of the multifunctionality of the ecological system in the target region within N years;

[0010] Step 3: Quantize the spatial pattern of the stability of the ecological system service supply and the spatial pattern of the ecological system resilience in the target region within N years, and extract the spatial distribution pattern of the stability of the ecological system in the target region within N years;

[0011] Step 4: According to the spatial distribution pattern of the multifunctionality of the ecological system and the spatial distribution pattern of the stability of the ecological system, extract the spatial distribution pattern of the key ecological zone of the ecological system in the target region within N years, and divide the key ecological zone into grades.

[0012] In step 2, the ecological system service includes soil conservation service, water production service, carbon fixation service and habitat maintenance service.

[0013] In step 2, the ecosystem multifunctionality index of the spatial distribution pattern of the target region in N years at the pixel scale is obtained by the following formula (I) and formula (II):

[0014] EMFI k = Mean(EMFI kj ) (I)

[0015]

[0016] wherein EMFI k represents the ecosystem multifunctionality index of the kth year;

[0017] EMFI kj represents the ecosystem multifunctionality index of the kth year at the threshold level;

[0018] F represents the total number of ecosystem services in the target region in N years;

[0019] i represents the ith ecosystem service;

[0020] f i represents the amount of the ith ecosystem service;

[0021] r i represents the conversion formula, when the f i of the ith ecosystem service is greater than the threshold t ij , the r i takes the value of 1, otherwise the r i takes the value of 0;

[0022] t ij represents the threshold of the ith ecosystem service at the threshold level;

[0023] Mean(x) represents the average value function.

[0024] In step 3, the ecosystem service supply stability index of the spatial pattern of the target region in N years at the pixel scale is obtained by the following formula (III):

[0025]

[0026] wherein ESSSI represents the ecosystem service supply stability index in N years;

[0027] f(x) represents the Jenks classification function;

[0028] n represents the number of ecosystem service supply in the target region in N years;

[0029] p represents the pth grid in the target region;

[0030] Mean pq represents the mean of the qth ecological service supply on the pth grid in N years;

[0031] SD pq represents the variance of the qth ecological service supply on the pth grid in N years.

[0032] wherein, in the step 3, the ecosystem resilience index of the spatial pattern of ecosystem resilience in the target region in N years at the pixel scale is obtained by the following formula (four):

[0033] ERI=F(Kendall τ(AR(1,NDVI),t)) (four)

[0034] wherein, ERI represents the ecosystem resilience index in N years;

[0035] f(x) represents the Jenks classification function;

[0036] Kendall τ represents the rank correlation coefficient, wherein one variable is AR(1,NDVI) and the other variable is time t;

[0037] AR(1,NDVI) represents the autoregressive coefficient of the N-year NDVI data set with a lag time of 1 month.

[0038] wherein, in the step 3, the ecosystem stability index of the spatial pattern of ecosystem stability in the target region in N years at the pixel scale is obtained by the following formula (five):

[0039] ESI=ω ESSSI ·f(ESSSI)+ω ERI ·f(ERI) (five)

[0040] wherein, ESI represents the ecosystem stability index in N years;

[0041] ESSSI represents the ecosystem service supply stability index in N years;

[0042] ω ESSSI represents the weight of ESSSI;

[0043] ERI represents the ecosystem resilience index in N years;

[0044] ω ERI represents the weight of ERI.

[0045] wherein, in the step 4, the spatial distribution pattern of the key ecological zone of the ecosystem in the target region in N years at the pixel scale is obtained by the following formula (six) and formula (seven):

[0046]

[0047] CEAs k = Top 30% (CEAI k ) (seven)

[0048] wherein CEAs k represents the kth year's ecosystem key ecological zone;

[0049] CEAI k represents the kth year's ecosystem key ecological zone identification index;

[0050] EMFI k represents the kth year's ecosystem multifunctionality index;

[0051] ESI represents the N-year ecosystem stability index;

[0052] represents a normalization function;

[0053] ω EMFI represents the weight of EMFI k ;

[0054] ω ESI represents the weight of ESI.

[0055] wherein in the step 4, in the key ecological zone grading process, the key ecological zone grade at the pixel scale is obtained by the following formula (eight):

[0056]

[0057] wherein L CEA represents the key ecological zone grade at each grid scale;

[0058] MEAN CEAI represents the CEAI average value at each grid scale in multiple years, i.e. the average value of the ecosystem key ecological zone identification index at each grid scale in multiple years;

[0059] SD CEAI represents the CEAI standard deviation at each grid scale in multiple years, i.e. the standard deviation of the ecosystem key ecological zone identification index at each grid scale in multiple years;

[0060] g(x) represents a natural break classification function.

[0061] wherein in the step 2, the j threshold value is respectively 25%, 50%, and 75%, obtaining three groups of the kth year's ecosystem multifunctionality index EMFI kThe average value is taken to further obtain the spatial distribution pattern of the ecosystem multifunctionality in the target region in N years.

[0062] The value of N is greater than or equal to 10, and the value of F is greater than or equal to 3.

[0063] The present application has the beneficial effects, including:

[0064] (1) The method for identifying the key ecological zone based on the ecosystem multifunctionality and stability provided by the present application can quickly and accurately identify the key ecological zone, which has important practical significance for maintaining the stability of regional ecosystem service supply, ensuring regional ecological safety, and realizing regional sustainable development.

[0065] (2) The method for identifying the key ecological zone based on the ecosystem multifunctionality and stability provided by the present application determines the pattern of the key ecological zone of the region from two angles of the ecosystem multifunctionality and stability, and the obtained key ecology is more accurate and reliable. BRIEF DESCRIPTION OF DRAWINGS

[0066] Figure 1 The overall logic diagram of the method for identifying the key ecological zone based on the ecosystem multifunctionality and stability of the present application is shown;

[0067] Figure 2a The spatial distribution pattern of the key ecological zone of the Yellow River Basin in 2005 obtained by the embodiment of the present application is shown;

[0068] Figure 2b The spatial distribution pattern of the key ecological zone of the Yellow River Basin in 2010 obtained by the embodiment of the present application is shown;

[0069] Figure 2c The spatial distribution pattern of the key ecological zone of the Yellow River Basin in 2015 obtained by the embodiment of the present application is shown;

[0070] Figure 2d The spatial distribution pattern of the key ecological zone of the Yellow River Basin in 2020 obtained by the embodiment of the present application is shown;

[0071] Figure 3 The schematic diagram of the classification of the key ecological zone of the Yellow River Basin obtained by the embodiment of the present application is shown;

[0072] Figure 4a The spatial distribution of the key ecological zone CEA and the existing important nature reserve ECA of the Yellow River Basin in 2005 in the embodiment of the present application is shown;

[0073] Figure 4b The spatial distribution of the key ecological zone CEA and the existing important nature reserve ECA of the Yellow River Basin in 2010 in the embodiment of the present application is shown;

[0074] Figure 4c Fig. 1 shows a schematic diagram of the spatial distribution of the CEA and ECA in the Yellow River Basin in 2015 according to an embodiment of the present application;

[0075] Figure 4d Fig. 2 shows a schematic diagram of the spatial distribution of the CEA and ECA in the Yellow River Basin in 2020 according to an embodiment of the present application. DETAILED DESCRIPTION

[0076] The application will be further described below by the accompanying drawings and embodiments. Through these descriptions, the features and advantages of the application will become more apparent.

[0077] The special word "exemplary" here means "as an example, embodiment or illustration". Any embodiment described as "exemplary" here is not necessarily interpreted as superior or better than other embodiments. Although various aspects of the embodiments are shown in the drawings, the drawings are not necessarily drawn to scale unless specifically indicated.

[0078] The application provides a key ecological zone identification method based on ecosystem multifunctionality and stability, as shown in Figure 1 The method comprises the following steps:

[0079] Step 1: Obtain a multi-source remote sensing data set in N years in the target area;

[0080] The multi-source remote sensing data set includes meteorological, soil, land use, terrain, leaf area index, net primary productivity and NDVI data;

[0081] The NDVI (Normalized Difference Vegetation Index) data can accurately reflect the surface vegetation coverage. In the present application, the NDVI data is obtained based on SPOT / VEGETATION and MODIS satellite remote sensing images. The NDVI data is applied in the research of vegetation dynamic change monitoring, land use / cover change detection, macroscopic vegetation cover classification and net primary productivity estimation in various scales.

[0082] In the present application, the actual meteorological data in N years in the target area is obtained by spatial interpolation through the daily scale data set of basic meteorological elements of China national ground meteorological stations, specifically including the spatio-temporal distribution pattern of air temperature, precipitation and wind speed.

[0083] In the present application, the soil data of the target area is obtained through the global public international soil reference and information center, specifically including the spatial distribution pattern of sand, silt, clay and organic carbon percentage in the soil.

[0084] In the present application, the spatio-temporal distribution pattern of land use in the target region for N years is obtained by the Institute of Geographic Science and Natural Resources Research, Chinese Academy of Sciences, and the land use categories mainly include forest land, grassland, water body, cultivated land, construction land and unused land.

[0085] In the present application, the spatial distribution pattern of elevation of the target region is obtained by the globally public ASTER DEM data set.

[0086] In the present application, the spatio-temporal distribution pattern of leaf area index, net primary productivity and NDVI of the target region for N years is obtained by the globally public MODIS data, which are MOD17A3 net primary productivity data set, MOD15A2H leaf area index data set and MOD13Q1 NDVI data set, respectively.

[0087] Step 2: quantifying the spatial distribution pattern of each ecosystem service in the target region for N years, and extracting the spatial distribution pattern of ecosystem multi-functionality in the target region for N years;

[0088] In step 2, the ecosystem services include soil conservation service, water production service, carbon fixation service and habitat maintenance service. In actual use, the selection of the type of ecological service can be adjusted according to the background status of the target region.

[0089] Preferably, the soil conservation service (SC) refers to the soil conservation amount per unit area per year, which reflects the ability to reduce soil erosion and land degradation, and represents the level of ecological function of protecting water and soil in the region. In the present application, the soil conservation amount is obtained at the pixel scale by the following formula (IX):

[0090] SC = ΔA = R x K x LS x (1 - C x P) (IX)

[0091] wherein ΔA represents the soil conservation amount (t / (hm 2 ·a));

[0092] R represents a rainfall erosion factor (MJ·mm / (hm 2 ·h·a);

[0093] K represents a soil erosion coefficient (t·hm 2 ·h / (MJ·mm·hm 2 );

[0094] LS represents a terrain factor (unitless);

[0095] C represents a vegetation coverage factor (unitless);

[0096] P represents a management factor (unitless);

[0097] In the present application, the specific values of R, K, LS, C and P are obtained by estimation or extraction from meteorological, soil, terrain, leaf area index, net primary productivity and NDVI data in step 1.

[0098] Habitat maintenance service (HQ) refers to the service provided by an ecosystem to maintain the life cycle of a species and genetic diversity, which belongs to the potential service to maintain the stability of an ecosystem. In the present application, the habitat quality is evaluated by using the modified habitat quality model to realize the spatial and temporal quantification of the habitat maintenance service of a region. On a pixel scale, the habitat maintenance service is obtained by the following formula (X) and formula (XI):

[0099]

[0100] wherein, MQ xp represents the modified habitat quality;

[0101] NDVI p represents the NDVI value in the pth grid;

[0102] D xp represents the stress degree of the pth grid in the land use type, i.e. the habitat degradation degree;

[0103] u represents a half-saturation constant;

[0104] z represents a model default parameter, preferably z = 2.5;

[0105] R represents the number of threat factors;

[0106] r represents a single threat factor;

[0107] Y r represents the number of grids occupied by the threat factor r;

[0108] y represents the grid where the threat factor r is located;

[0109] W r represents the weight of the threat factor r;

[0110] r y represents the threat factor intensity;

[0111] i rxy represents the influence of the threat factor on each grid;

[0112] S pr represents the relative sensitivity of the land use type to different threat factors.

[0113] Water yield service (WY) refers to the water yield per unit area of land cover in a certain period, which specifically includes surface runoff, interflow, canopy interception and litter water retention. Water yield is closely related to precipitation, soil, topography, surface cover type, root depth, soil depth and plant transpiration, etc. At the pixel scale, it is obtained by the following formula (twelve):

[0114]

[0115] wherein Y hp represents the annual water yield (m 3 ) of the pth grid on the hth land use type;

[0116] AET ph represents the annual average evapotranspiration (mm) of the pth grid on the hth land use type;

[0117] P p represents the annual average rainfall (mm) of the pth grid;

[0118] is the ratio coefficient of actual evaporation and average annual precipitation, which is obtained by using the approximate algorithm based on Budyko curve.

[0119] Carbon sequestration service (CS) refers to the ability of vegetation to absorb carbon dioxide in the atmosphere. As an important regulating service, it plays an important role in slowing the rise of atmospheric CO2 concentration, regulating global climate and maintaining global carbon balance. At the pixel scale, it is obtained by the following formula (thirteen):

[0120] CS=C total =(C above +C below +C soil +C dead )×S (thirteen)

[0121] wherein C total represents the total carbon storage on the estimation unit;

[0122] C above represents the carbon density in aboveground biomass;

[0123] C below represents the carbon density in belowground biomass;

[0124] C soil represents the carbon density in soil;

[0125] C dead represents the carbon density in dead organic matter;

[0126] S represents the area of the estimation unit.

[0127] Preferably, in step 2, the ecosystem multifunctionality index of the spatial distribution pattern of ecosystem multifunctionality in the target region in N years at pixel scale is obtained by the following formula (I) and formula (II):

[0128] EMFI k = Mean(EMFI kj ) (I)

[0129]

[0130] wherein EMFI k represents the ecosystem multifunctionality index of the kth year;

[0131] EMFI kj represents the ecosystem multifunctionality index of the kth year at the threshold level of i;

[0132] F represents the total number of ecosystem services in the target region in N years;

[0133] i represents the ith ecosystem service;

[0134] f i represents the amount of the ith ecosystem service;

[0135] r i represents the conversion formula, when the f i of the ith ecosystem service is greater than the threshold t ij , the r i takes the value of 1, otherwise the r i takes the value of 0;

[0136] t ij represents the threshold of the ith ecosystem service at the jth threshold level, which is obtained by multiplying the maximum value by the corresponding threshold;

[0137] Mean(x) represents the average value function.

[0138] Preferably, in step 2, the jth threshold takes the value of 25%, 50%, and 75% respectively, to obtain three groups of the ecosystem multifunctionality index EMFI k of the kth year, and the average value thereof is used to further obtain the spatial distribution pattern of ecosystem multifunctionality in the target region in N years.

[0139] In a preferred embodiment, the value of N is greater than or equal to 10, and the value of F is greater than or equal to 3.

[0140] In the present application, the specific value of N is not specifically limited, and those skilled in the art can determine it according to the location of the target area and other conditions. If the target area has a large activity intensity and a large human activity disturbance, the value of N can be appropriately increased to improve the reliability of the calculation result. If the target area has a small activity intensity and a small human activity disturbance, the value of N can be appropriately reduced to reduce the calculation workload.

[0141] Step 3: Quantifying the spatial pattern of the stability of ecosystem service supply and the spatial pattern of the ecosystem resilience in the target area in N years, and extracting the spatial distribution pattern of the stability of the ecosystem in the target area in N years;

[0142] In the step 3, the spatial pattern of the stability of ecosystem service supply in the target area in N years is the ecosystem service supply stability index at the pixel scale, which is obtained by the following formula (three):

[0143]

[0144] Wherein, ESSSI represents the ecosystem service supply stability index in N years, and the larger the value is, the more stable the service supply is;

[0145] f(x) represents the Jenks classification function;

[0146] n represents the number of ecosystem service supply in the target area in N years;

[0147] p represents the pth grid in the target area; the grid refers to a spatial data model defined by a series of equal-sized pixels, in which the data is arranged by rows and columns of pixels, and each pixel contains an attribute value and a set of position coordinates;

[0148] Mean pq represents the mean of the qth ecosystem service supply in the pth grid in N years;

[0149] SD pq represents the variance of the qth ecosystem service supply in the pth grid in N years.

[0150] Preferably, in the step 3, the spatial pattern of the ecosystem resilience in the target area in N years is the ecosystem resilience index at the pixel scale, which is obtained by the following formula (four):

[0151] ERI=f(Kendallτ(AR(1,NDVI),t)) (four)

[0152] Wherein, ERI represents the ecosystem resilience index in N years, and the larger the value is, the stronger the ecosystem resilience is;

[0153] f(x) represents the Jenks classification function;

[0154] Kendall τ represents the rank correlation coefficient, where one variable is AR(1, NDVI) and the other variable is time t;

[0155] AR(1, NDVI) represents the autoregressive coefficient of the N-year NDVI dataset with a lag time of 1 month.

[0156] Preferably, in the step 3, the spatial distribution pattern of the ecosystem stability in the target region in N years is an ecosystem stability index at the pixel scale, which is obtained by the following formula (V):

[0157] ESI = ω ESSSI · f(ESSSI) + ω ERI · f(ERI) (V)

[0158] Wherein, ESI represents the ecosystem stability index in N years;

[0159] ESSSI represents the ecosystem service supply stability index in N years;

[0160] ERI represents the ecosystem resilience index in N years;

[0161] ω ESSSI represents the weight of ESSSI;

[0162] ω ERI represents the weight of ERI.

[0163] In the present application, the entropy weight method is used to calculate the weight of different indexes. Preferably, in order to reduce subjectivity and balance the importance of each index, the value of ω ESSSI may be 0.5, and the value of ω ERI may be 0.5.

[0164] Step 4: According to the spatial distribution pattern of the ecosystem multi-functionality and the spatial distribution pattern of the ecosystem stability, the spatial distribution pattern of the key ecological zone of the ecosystem in the target region in N years is extracted, and the key ecological zone is classified.

[0165] In the step 4, the spatial distribution pattern of the key ecological zone of the ecosystem in the target region in N years is the key ecological zone of the ecosystem at the pixel scale, which is obtained by the following formula (VI) and formula (VII):

[0166]

[0167] CEAs k = Top 30% (CEAI k ) (VII)

[0168] Wherein, CEAs kEESI k represents the ecosystem key ecological zone in the kth year;

[0169] CEAI k EESI k represents the ecosystem key ecological zone identification index in the kth year; that is, the area corresponding to the top 30% of the largest EESI k in the kth year is the key ecological zone.

[0170] EMFI k EMFI k represents the ecosystem multifunctionality index in the kth year;

[0171] ESI represents the N-endophyte ecosystem stability index;

[0172] f represents a normalization function;

[0173] ω EMFI ω k represents the weight of EMFI k; k

[0174] ω ESI ω represents the weight of ESI.

[0175] In the present application, the entropy weight method is used to calculate the weight of different indexes. Preferably, in order to reduce subjectivity and balance the importance of each index, the value of ω k is 0.5, and the value of ω is 0.5. EMFI ESI

[0176] In a preferred embodiment, in step 4, in the process of classifying the key ecological zones, the key ecological zone level is obtained at the pixel scale by the following formula (eight):

[0177]

[0178] wherein L k represents the key ecological zone level at each grid scale; CEA

[0179] MEAN CEAI MEAN k represents the average value of CEAI in multiple years at each grid scale, i.e., the average value of the ecosystem key ecological zone identification index in multiple years at each grid scale;

[0180] SD CEAI SD k represents the standard deviation of CEAI in multiple years at each grid scale, i.e., the standard deviation of the ecosystem key ecological zone identification index in multiple years at each grid scale;

[0181] g(x) represents a natural break classification function.

[0182] Embodiment

[0183] The Yellow River Basin is selected as the target region, and the key ecological zone of the target region is identified. The specific identification process includes the following steps:​​​​

[0184] Step 1: Obtain multi-source remote sensing data sets of the target region in N years;

[0185] The multi-source remote sensing data sets include meteorological, soil, land use, terrain, leaf area index, net primary productivity, and NDVI data; N is 15 years;

[0186] Among them, each multi-source remote sensing data source is shown in Table 1 below;

[0187] Table 1 Data Description

[0188]

[0189]

[0190] Step 2: Quantify the spatial distribution pattern of each ecosystem service in the target region in N years, and extract the spatial distribution pattern of ecosystem multi-functionality in the target region in N years;

[0191] The ecosystem services include soil conservation services, water production services, carbon sequestration services, and habitat maintenance services.

[0192] In step 2, the spatial distribution pattern of ecosystem multi-functionality in the target region in N years is the ecosystem multi-functionality index at the pixel scale, which is obtained by the following formula (I) and formula (II):

[0193] EMFI k = Mean(EMFI kj ) (I)

[0194]

[0195] Among them, EMFI k represents the ecosystem multi-functionality index of the kth year;

[0196] EMFI kj represents the ecosystem multi-functionality index of the kth year at the j threshold level;

[0197] F represents the total number of ecosystem services in the target region in N years; the specific value is 4;

[0198] i represents the ith ecosystem service;

[0199] f i represents the amount of the ith ecosystem service;

[0200] r i represents the conversion formula, when the f i of the ith ecosystem service is greater than the threshold t ij , the r iis 1, otherwise r i is 0.

[0201] t ij denotes the threshold value of the ith ecosystem service at the threshold level;

[0202] Mean(x) represents the mean function.

[0203] In step 2, the j-threshold values are 25%, 50%, and 75%, respectively, to obtain three groups of ecosystem multifunctionality index EMFI in the k-th year k , and the average value is taken to further obtain the spatial distribution pattern of ecosystem multifunctionality in the target area N years.

[0204] Step 3: Quantify the spatial pattern of ecosystem service supply stability and ecosystem resilience in the target area N years, and extract the spatial distribution pattern of ecosystem stability in the target area N years;

[0205] In the step 3, the spatial pattern of ecosystem service supply stability in the target area N years is the ecosystem service supply stability index at the pixel scale, which is obtained by the following formula (three):

[0206]

[0207] Wherein, ESSSI represents the ecosystem service supply stability index in N years, and the larger the value represents the more stable the service supply;

[0208] f(x) represents the Jenks classification function;

[0209] n represents the number of ecosystem service supply in the target area N years;

[0210] p represents the pth grid in the target area;

[0211] Mean pq represents the mean value of the qth ecosystem service supply in the pth grid in N years;

[0212] SD pq represents the variance of the qth ecosystem service supply in the pth grid in N years.

[0213] In the step 3, the spatial pattern of ecosystem resilience in the target area N years is the ecosystem resilience index at the pixel scale, which is obtained by the following formula (four):

[0214] ERI=F(Kendall tau(AR(1, NDVI), t)) (four)

[0215] wherein ERI represents the ecosystem resilience index in N years, and the greater the value, the stronger the ecosystem resilience;

[0216] f(x) represents the Jenks classification function;

[0217] Kendall τ represents the rank correlation coefficient, wherein one variable is AR(1, NDVI) and the other variable is time t;

[0218] AR(1, NDVI) represents the autoregressive coefficient of the NDVI data set in N years with a lag time of 1 month.

[0219] In the step 3, the spatial distribution pattern of the ecosystem stability in the target region in N years is the ecosystem stability index at the pixel scale, which is obtained by the following formula (V):

[0220] ESI = ω ESSSI · F(ESSSI) + ω ERI · f(ERI) (V)

[0221] wherein ESI represents the ecosystem stability index in N years;

[0222] ESSSI represents the ecosystem service supply stability index in N years;

[0223] ERI represents the ecosystem resilience index in N years;

[0224] ω ESSSI represents the weight of ESSSI, and the value is 0.5;

[0225] ω ERI represents the weight of ERI, and the value is 0.5.

[0226] Step 4: According to the spatial distribution pattern of the ecosystem multi-functionality and the spatial distribution pattern of the ecosystem stability, the spatial distribution pattern of the key ecological zone of the ecosystem in the target region in N years is extracted, and the key ecological zone is classified.

[0227] In the step 4, the spatial distribution pattern of the key ecological zone of the ecosystem in the target region in N years is the key ecological zone of the ecosystem at the pixel scale, which is obtained by the following formula (VI) and formula (VII):

[0228]

[0229] CEAs k = Top30% (CEAI k ) (VII)

[0230] wherein CEAs k represents the key ecological zone of the ecosystem in the kth year;

[0231] CEAI k denotes the ecosystem key ecological zone identification index in the kth year;

[0232] EMFI k denotes the ecosystem multi-functionality index in the kth year;

[0233] ESI denotes the ecosystem stability index in N years;

[0234] denotes a normalization function;

[0235] ω EMFI denotes the weight of EMFI k , and is 0.5;

[0236] ω ESI denotes the weight of ESI, and is 0.5.

[0237] In the step 4, in the key ecological zone grading process, the key ecological zone grade is obtained at the pixel scale by the following formula (eight):

[0238]

[0239] wherein L CEA denotes the key ecological zone grade at each grid scale;

[0240] MEAN CEAI denotes the average value of CEAI in multiple years at each grid scale;

[0241] SD CEAI denotes the standard deviation of CEAI in multiple years at each grid scale;

[0242] g(x) denotes a natural break classification function.

[0243] The finally obtained key ecological zone distribution and grading of the Yellow River Basin in 2005-2020 are shown in FIG. 2 and FIG. 3. Figure 3 Among them, Figures 2a to 2d shows the spatial distribution pattern of the key ecological zones of the Yellow River Basin in 2005-2020; Figure 3 shows the key ecological zone grading of the Yellow River Basin.

[0244] Further, the above identification result is superimposed and compared with the most important nature reserves ECAs provided by the Ministry of Ecology and Environment, specifically, the range of ECAs includes the priority areas for biodiversity conservation PABCs and the national nature reserves NNRs.

[0245] The above example result shows that among the existing 28 ECAs, 23 are located within the range of the CEA identified in the above example, such asFigures 4a to 4d The recognition accuracy reaches 82%.

[0246] Except for the West Ordos-Helanshan-Yinshan OHYZ, most of the most important nature reserves in other regions are recognized as CEAs by the method in the above embodiment. The main reason for the low recognition accuracy of the West Ordos-Helanshan-Yinshan OHYZ is that the OHYZ mainly focuses on the protection of desert ecosystems and geological relics, such as the E'tuoqi Dinosaur Fossil National Nature Reserve and the Inner Mongolia West Ordos National Nature Reserve. Since the services and functions of the desert ecosystem are not considered in the above embodiment, some of the most important nature reserves in the OHYZ are not recognized as CEAs by the above embodiment.

[0247] In summary, the results obtained by using the key ecological zone recognition method based on the ecosystem multifunctionality and stability proposed in the present patent are reliable, and have an important guiding role for determining the priority protection areas in the Yellow River Basin. The recognition framework can realize the rapid recognition of key ecological zones in a large-scale range.

[0248] The above describes the present application in combination with preferred embodiments, but these embodiments are only exemplary and serve only to illustrate. On this basis, various substitutions and improvements can be made to the present application, and these all fall within the protection scope of the present application.

Claims

1. A method for identifying key ecological regions based on ecosystem multifunctionality and stability, characterized by: The method comprises the following steps: Step 1: Obtain a multi-source remote sensing dataset of the target area within N years; Step 2: Quantify the spatial distribution pattern of each ecosystem service in the target area within N years and extract the spatial distribution pattern of ecosystem multifunctionality in the target area within N years; Step 3: Quantify the spatial pattern of ecosystem service supply stability and ecosystem resilience in the target area over N years, and extract the spatial distribution pattern of ecosystem stability in the target area over N years; Step 4: Based on the spatial distribution pattern of ecosystem multifunctionality and the spatial distribution pattern of ecosystem stability, extract the spatial distribution pattern of key ecological zones of the target area within N years and classify the key ecological zones; In step 2, the ecosystem multifunctionality index of the target area's spatial distribution pattern of ecosystem multifunctionality at the pixel scale over N years is obtained using the following formulas (1) and (2): EMFI k = Mean(EMFI kj ) (1) Among them, EMFI k represents the ecosystem multifunctionality index in year k; EMFI kj represents the ecosystem multifunctionality index at the jth threshold level in year k; F represents the total number of ecosystem services in the target area within N years; i represents the i-th ecosystem service; f i represents the amount of the i-th ecosystem service; r i Represents the conversion formula, when the f of the i-th ecosystem service i is greater than the threshold t ij When the r i The value is 1, otherwise r i The value is 0; t ij represents the threshold value of the i-th ecosystem service at the j-th threshold level; Mean(x) represents the mean value function; In step 3, the ecosystem stability index of the spatial distribution pattern of ecosystem stability in the target area within N years at the pixel scale is obtained by the following formula (5): ESI = ω ESSSI ·f(ESSSI) + ω ERI ·f(ERI) (Five) Among them, ESI represents the ecosystem stability index within N years; ESSSI represents the ecosystem service supply stability index within N years; ω ESSSI represents the weight of ESSSI; ERI represents the ecosystem resilience index within N years; ω ERI represents the weight of ERI; In step 4, the spatial distribution pattern of the key ecological zones of the ecosystem in the target area within N years is obtained at the pixel scale by the following formulas (6) and (7): CEAs k = Top30% (CEAI k ) (Seven) Among them, CEAs k represents the key ecological area of ​​the ecosystem in year k; CEAI k represents the key ecological zone identification index of the ecosystem in year k; EMFI k represents the ecosystem multifunctionality index in year k; ESI represents the ecosystem stability index within N years; represents the normalization function; ω EMFI Indicates EMFI k The weight of ω ESI Indicates the weight of ESI.

2. The method for identifying key ecological areas based on ecosystem multifunctionality and stability according to claim 1, characterized in that: In step 2, the ecosystem services include soil conservation services, water production services, carbon sequestration services and habitat maintenance services.

3. The method for identifying key ecological areas based on ecosystem multifunctionality and stability according to claim 1, characterized in that: In step 3, the ecosystem service supply stability index of the spatial pattern of ecosystem service supply stability in the target area over N years at the pixel scale is obtained by the following formula (3): Among them, ESSSI represents the ecosystem service supply stability index within N years; f(x) represents the Jenks classification function; n represents the quantity of ecosystem services provided in the target area within N years; p represents the pth grid in the target area; Mean pq represents the mean value of the qth ecological service supply on the pth grid over N years; SD pq It represents the variance of the qth ecological service supply on the pth grid over N years.

4. The method for identifying key ecological areas based on ecosystem multifunctionality and stability according to claim 1, characterized in that: In step 3, the ecosystem resilience index of the target area's ecosystem resilience spatial pattern at the pixel scale over N years is obtained using the following formula (IV): ERI=f(Kendallτ(AR(1,NDVI),t)) (4) Wherein, ERI represents the ecosystem resilience index within N years; f(x) represents the Jenks classification function; Kendallτ represents the rank correlation coefficient, where one variable is AR(1, NDVI) and the other variable is time t; AR(1,NDVI) represents the autoregressive coefficient of the N-year NDVI dataset with a lag time of 1 month.

5. The method for identifying key ecological areas based on ecosystem multifunctionality and stability according to claim 1, characterized in that: In step 4, during the classification of key ecological zones, the key ecological zone levels are obtained at the pixel scale using the following formula (VIII): Among them, L CEA Indicates the level of key ecological areas at each grid scale; MEAN CEAI It represents the average value of CEAI over many years at each grid scale, that is, the average value of the key ecological zone identification index of ecosystems over many years at each grid scale; SD CEAI It represents the standard deviation of CEAI over multiple years at each grid scale, that is, the standard deviation of the key ecological zone identification index of ecosystems over multiple years at each grid scale; g(x) represents the natural fracture classification function.

6. The method for identifying key ecological areas based on ecosystem multifunctionality and stability according to claim 1, characterized in that: In step 2, the j threshold is set to 25%, 50%, and 75%, respectively, and three sets of ecosystem multifunctionality indices EMFI for the kth year are obtained: k , and take its average value to further obtain the spatial distribution pattern of ecosystem multifunctionality in the target area within N years.

7. The method for identifying key ecological areas based on ecosystem multifunctionality and stability according to claim 1, characterized in that: The value of N is greater than or equal to 10, and the value of F is greater than or equal to 3.

Images