Key ecological region identification method based on multifunctionality and stability of ecological system

Abstract

The invention discloses a key ecological region identification method based on the multifunctionality and stability of an ecological system, and the method combines the multifunctionality and stability of the ecological system to determine the key ecological region pattern of a region, and enables limited resources and funds to be input into the key ecological region of the region. Therefore, the maximization of regional ecological protection benefits is realized.

Description

Technical Field

[0001] The invention relates to the technical fields of regional ecosystem governance and regional land use planning and management, and in particular to a key ecological zone identification method based on the multifunctionality and stability of the ecosystem. Background Art

[0002] Critical Ecological Areas (CEAs) are important areas for maintaining biodiversity and ecosystem functions, continuously providing ecosystem services to human society, and playing an important role in regional ecological security. However, in recent years, under the dual impact of climate change and human activities, the functional degradation of regional critical ecological areas has become increasingly severe, and the problem of declining ecosystem service supply capacity cannot be ignored. In the process of regional ecological protection policy formulation and restoration, due to the limitation of investment funds and resources, determining the order of protection, identifying regional critical ecological areas and priority protection hotspots are the most basic and critical stages. Investing limited resources and funds in regional critical ecological areas is conducive to improving the benefits of ecological protection and management. Therefore, the identification of regional critical ecological areas has become an important tool for improving ecosystem degradation and maintaining reasonable landscape planning. At present, related research has received widespread attention worldwide. Therefore, in-depth research on the identification of regional critical ecological areas and the formulation of reasonable regional spatial planning plans for regional planning and development will provide a fundamental solution for coordinating regional ecological protection and sustainable social and economic development.

[0003] Ecosystem services are the bridge between natural systems and human well-being. Ecosystems continuously provide products and services to humans through ecological functions to meet human well-being needs. At present, the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services (IPBES) pointed out in its latest survey report that the ability of nature to maintain its contribution to a 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 zones based on ecosystem services. Key ecological zones based on ecosystem services are defined as "areas with important ecological service function supply compared to other regions", which should receive full attention and protection. From the perspective of ecosystem services, the identification of regional key ecological zones has received increasing attention, 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 the static assessment of factors is mainly based on, such as the spatiotemporal superposition analysis of multiple single ecological service functions based on 3S technology and model simulation, focusing on the static assessment of ecosystem service supply or long-term mean assessment. Such methods pay too much attention to the static or mean characteristics of ecosystems, and insufficient consideration of the dynamic change characteristics of long-term ecosystems. As the dynamic stability and multifunctionality of long-term ecosystems are important intrinsic characteristics of ecosystems, they are conducive to the stable and efficient supply of regional ecological services. However, related studies have not systematically introduced the concepts and determination methods related to ecosystem multifunctionality and stability, and their inclusion in the framework of regional key ecological zone identification needs to be further deepened.

[0005] Based on the above problems, the inventors have conducted in-depth research on the dynamic stability and multifunctionality of long-term ecosystems, hoping to discover a method to quickly and accurately obtain key ecological zones based on this. Summary of the invention

[0006] In order to overcome the above problems, the inventors conducted intensive research and designed a method for identifying key ecological zones based on the multifunctionality and stability of ecosystems. In this method, the key ecological zone pattern of a region is determined from the two perspectives of ecosystem multifunctionality and stability, and limited resources and funds are invested in regional key ecological zones, thereby maximizing the benefits of regional ecological protection, thereby completing the present invention.

[0007] Specifically, the present invention aims to provide a method for identifying key ecological areas based on the multifunctionality and stability of ecosystems, which method comprises the following steps:

[0008] Step 1: Obtain multi-source remote sensing datasets of the target area within N years;

[0009] 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;

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

[0011] 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.

[0012] Wherein, in step 2, the ecosystem services include soil conservation services, water production services, carbon sequestration services and habitat maintenance services.

[0013] Among them, in step 2, the ecosystem multifunctionality index of the spatial distribution pattern of ecosystem multifunctionality in the target area within N years at the pixel scale is obtained by the following formulas (I) and (II):

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

[0015]

[0016] Among them, EMFI k represents the ecosystem multifunctionality index in the kth year;

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

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

[0019] i represents the i-th ecosystem service;

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

[0021] 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;

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

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

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

[0025]

[0026] Among them, ESSSI represents the ecosystem service supply stability index within N hands;

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

[0028] n represents the quantity of ecosystem service supply in the target area within N years;

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

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

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

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

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

[0034] Among them, ERI represents the ecosystem resilience index within N hands;

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

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

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

[0038] Among them, in the 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 (V):

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

[0040] Among them, ESI represents the ecosystem stability index within N years;

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

[0042] ω ESSSI represents the weight of ESSSI;

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

[0044] ω ERI Indicates the weight of ERI.

[0045] Among them, in the step 4, the spatial distribution pattern of the key ecological zones of the ecosystem in the target area within N years is obtained by the following formula (VI) and formula (VII) at the pixel scale:

[0046]

[0047] CEAs k =Top30%(CEAI k ) (seven)

[0048] Among them, CEAs k represents the key ecological area of ​​the ecosystem in year k;

[0049] CEAI k represents the key ecological zone identification index of the ecosystem in year k;

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

[0051] ESI represents the ecosystem stability index within N years;

[0052] represents the normalization function;

[0053] ω EMFI Represents EMFI k The weight of

[0054] ω ESI Indicates the weight of ESI.

[0055] Wherein, in the step 4, during the classification of the key ecological areas, at the pixel scale, the key ecological area level is obtained by the following formula (VIII):

[0056]

[0057] Among them, L CEA Indicates the level of key ecological areas at each grid scale;

[0058] 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;

[0059] 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;

[0060] g(x) represents the natural fracture classification function.

[0061] In step 2, the j threshold is set to 25%, 50%, and 75%, respectively, and the three sets of ecosystem multifunctionality index EMFI for the kth year are obtained: k, and take the average value to further obtain the spatial distribution pattern of ecosystem multifunctionality in the target area within N years.

[0062] Among them, 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 beneficial effects of the present invention include:

[0064] (1) The key ecological zone identification method based on ecosystem multifunctionality and stability provided by the present invention can quickly and accurately identify key ecological zones, which has important practical significance for maintaining the stability of regional ecosystem service supply, ensuring regional ecological security, and achieving regional sustainable development;

[0065] (2) According to the key ecological zone identification method based on ecosystem multifunctionality and stability provided by the present invention, the method determines the key ecological zone pattern of a region by combining the two perspectives of ecosystem multifunctionality and stability, and the key ecology obtained is more accurate and reliable. BRIEF DESCRIPTION OF THE DRAWINGS

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

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

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

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

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

[0071] Figure 3 A schematic diagram showing the classification of key ecological zones in the Yellow River Basin obtained in an embodiment of the present application is shown;

[0072] Figure 4a The figure shows the spatial distribution diagram of the key ecological zones CEAs and the existing important nature reserves ECAs in the Yellow River Basin in 2005 in the embodiment of the present application;

[0073] Figure 4b The figure shows the spatial distribution diagram of the key ecological zones CEAs and the existing important nature reserves ECAs in the Yellow River Basin in 2010 in the embodiment of the present application;

[0074] Figure 4c The figure shows the spatial distribution diagram of the key ecological zones CEAs and the existing important nature reserves ECAs in the Yellow River Basin in 2015 in the embodiment of the present application;

[0075] Figure 4d A schematic diagram of the spatial distribution of key ecological zones (CEAs) and existing important nature reserves (ECAs) in the Yellow River Basin in 2020 in an embodiment of the present application is shown. DETAILED DESCRIPTION

[0076] The present invention will be further described in detail below through the accompanying drawings and embodiments. Through these descriptions, the characteristics and advantages of the present invention will become more clear and distinct.

[0077] The word "exemplary" is used exclusively herein to mean "serving as an example, embodiment, or illustration." Any embodiment described herein as "exemplary" is not necessarily to be construed as preferred or advantageous over other embodiments. Although various aspects of the embodiments are shown in the drawings, the drawings are not necessarily drawn to scale unless otherwise noted.

[0078] The present invention provides a method for identifying key ecological areas based on the multifunctionality and stability of ecosystems, such as Figure 1 As shown in , the method comprises the following steps:

[0079] Step 1: Obtain multi-source remote sensing datasets of the target area within N years;

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

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

[0082] In this application, the daily-scale dataset of basic meteorological elements of China's national ground meteorological stations is used to obtain the actual meteorological data of the target area within N years through spatial interpolation, including the spatiotemporal distribution patterns of temperature, precipitation, and wind speed.

[0083] In this application, soil data of the target area are obtained through the globally open 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 this application, the spatiotemporal distribution pattern of land use in the target area over N years was obtained through the Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences. The land use categories mainly include: forest land, grassland, water body, cultivated land, construction land and unused land.

[0085] In this application, the spatial distribution pattern of the elevation of the target area is obtained through the globally public ASTER DEM dataset.

[0086] In this application, the spatiotemporal distribution patterns of leaf area index, net primary productivity, and NDVI in the target area for N years were obtained through the globally public MODIS data, which are the MOD17A3 net primary productivity dataset, the MOD15A2H leaf area index dataset, and the MOD13Q1 NDVI dataset.

[0087] 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;

[0088] In step 2, the ecosystem services include soil conservation services, water production services, carbon sequestration services and habitat maintenance services. In actual use, the selection of ecosystem service types can be adjusted according to the background conditions of the target area.

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

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

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

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

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

[0094] LS represents the topographic factor (unitless);

[0095] C represents the vegetation cover factor (unitless);

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

[0097] In this application, the specific values ​​of R, K, LS, C and P are estimated or extracted by the meteorological, soil, topographic, leaf area index, net primary productivity and NDVI data in step 1.

[0098] Habitat maintenance service (HQ) refers to the services provided by ecosystems to maintain species life cycles and genetic diversity, and is a potential service for maintaining ecosystem stability. In this application, the modified habitat quality model is used to evaluate habitat quality and achieve spatiotemporal quantification of regional habitat maintenance services, which is obtained at the pixel scale through the following formulas (10) and (11):

[0099]

[0100] Among them, MQ xp represents the corrected habitat quality;

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

[0102] D xp It represents the degree of stress of the pth grid in the land use type, that is, the degree of habitat degradation;

[0103] u represents the half-saturation constant;

[0104] z represents the default parameters of the model, preferably z=2.5;

[0105] R represents the number of threat factors;

[0106] r represents a single threat factor;

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

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

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

[0110] r y Indicates the intensity of threat factors;

[0111] i rxy Indicates the impact of threat factors on each grid;

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

[0113] Water production service (WY) refers to the water production per unit area of ​​ground cover in a certain period of time, including surface runoff, soil flow, canopy interception and litter water retention. Water production is closely related to precipitation, soil, topography, ground cover type, root depth, soil depth and plant transpiration, and can be obtained at the pixel scale by the following formula (XII):

[0114]

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

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

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

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

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

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

[0121] Among them, C total represents the total carbon stock 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 Indicates the carbon density in the 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 area within N years at the pixel scale is obtained by the following formulas (I) and (II):

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

[0129]

[0130] Among them, EMFI k represents the ecosystem multifunctionality index in the kth year;

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

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

[0133] i represents the i-th ecosystem service;

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

[0135] r i represents the conversion formula, when the f of the i-th ecosystem service i is greater than the threshold t ij When r i The value is 1, otherwise r i The value is 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 mean function.

[0138] Preferably, in step 2, the j threshold is set to 25%, 50%, and 75%, respectively, to obtain three sets of ecosystem multifunctionality index EMFI for the kth year: k , and take the average value to further obtain the spatial distribution pattern of ecosystem multifunctionality in the target area within 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 this 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, etc. If the human activity intensity in the target area is large and the human activity disturbance is large, the N value can be appropriately increased to improve the reliability of the calculation result. If the human activity intensity in the target area is small and the human activity disturbance is small, the N value can be appropriately reduced to reduce the calculation workload.

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

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

[0143]

[0144] Among them, ESSSI represents the ecosystem service supply stability index within N hands, and the larger the value, the more stable the service supply;

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

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

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

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

[0149] SD pq It represents the variance of the qth ecological service supply on the pth grid within N years.

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

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

[0152] Among them, ERI represents the ecosystem resilience index within N years, and the larger the value, the stronger the ecosystem resilience;

[0153] f(x) represents 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 autoregression coefficient of the N-year NDVI dataset with a lag time of 1 month.

[0156] Preferably, 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 (V):

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

[0158] Among them, ESI represents the ecosystem stability index within N years;

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

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

[0161] ω ESSSI represents the weight of ESSSI;

[0162] ω ERI Indicates the weight of ERI.

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

[0164] 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.

[0165] In step 4, the spatial distribution pattern of the key ecological zones of the ecosystem within the target area N hands is obtained at the pixel scale by the following formula (VI) and formula (VII):

[0166]

[0167] CEAs k =Top30%(CEAI k )(seven)

[0168] Among them, CEAs krepresents the key ecological area of ​​the ecosystem in year k;

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

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

[0171] ESI stands for ecosystem stability index within N hands;

[0172] represents the normalization function;

[0173] ω EMFI Represents EMFI k The weight of

[0174] ω ESI Indicates the weight of ESI.

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

[0176] In a preferred embodiment, in step 4, during the classification of key ecological areas, the key ecological area level is obtained at the pixel scale by the following formula (VIII):

[0177]

[0178] Among them, L CEA Indicates the level of key ecological areas at each grid scale;

[0179] 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;

[0180] 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;

[0181] g(x) represents the natural fracture classification function.

[0182] Example

[0183] The Yellow River Basin was selected as the target area for the identification of key ecological zones. The specific identification process included the following steps:

[0184] Step 1: Obtain multi-source remote sensing datasets of the target area within N years;

[0185] The multi-source remote sensing data set includes meteorological, soil, land use, topography, leaf area index, net primary productivity and NDVI data; the value of N is 15 years;

[0186] Among them, the data sources of various multi-source remote sensing data are 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 area within N years, and extract the spatial distribution pattern of ecosystem multifunctionality in the target area within 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 ecosystem multifunctionality index of the target area in the pixel scale for the spatial distribution pattern of ecosystem multifunctionality within N years is obtained by the following formulas (I) and (II):

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

[0194]

[0195] Among them, EMFI k represents the ecosystem multifunctionality index in the kth year;

[0196] EMFI kj represents the ecosystem multifunctionality index at the jth threshold level in the kth year;

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

[0198] i represents the i-th ecosystem service;

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

[0200] 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 iThe value is 1, otherwise r i The value is 0;

[0201] t ij represents the threshold of the ith ecosystem service under the threshold level;

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

[0203] In step 2, the j threshold is set to 25%, 50%, and 75%, respectively, and three sets of ecosystem multifunctionality index EMFI for the kth year are obtained: k , and take the average value to further obtain the spatial distribution pattern of ecosystem multifunctionality in the target area within N years.

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

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

[0206]

[0207] Among them, ESSSI represents the ecosystem service supply stability index within N hands, and the larger the value, the more stable the service supply;

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

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

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

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

[0212] SD pq It represents the variance of the qth ecological service supply on the pth grid within N years.

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

[0214] ERI=F(Kendallτ(AR(1,NDVI),t)) (4)

[0215] Among them, ERI represents the ecosystem resilience index within N years, and the larger the value, the stronger the ecosystem resilience;

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

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

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

[0219] 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 (V):

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

[0221] Among them, ESI represents the ecosystem stability index within N hands;

[0222] ESSSI represents the ecosystem service supply stability index within N hands;

[0223] ERI stands for Ecosystem Resilience Index within N hands;

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

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

[0226] 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.

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

[0228]

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

[0230] Among them, CEAs k represents the key ecological area of ​​the ecosystem in year k;

[0231] CEAI k represents the key ecological zone identification index of the ecosystem in year k;

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

[0233] ESI represents the ecosystem stability index within N years;

[0234] represents the normalization function;

[0235] ω EMFI Represents EMFI k The weight of is 0.5;

[0236] ω ESI Indicates the weight of ESI, and its value is 0.5.

[0237] In step 4, during the classification of key ecological areas, the key ecological area level is obtained at the pixel scale by the following formula (VIII):

[0238]

[0239] Among them, L CEA Indicates the level of key ecological areas at each grid scale;

[0240] MEAN CEAI It represents the average CEAI value over multiple years at each grid scale;

[0241] SD CEAI It represents the standard deviation of CEAI between multiple years at each grid scale;

[0242] g(x) represents the natural fracture classification function.

[0243] The final distribution and classification of key ecological zones in the Yellow River Basin from 2005 to 2020 are shown in Figures 2 and Figure 3 As shown. Among them, Figure 2a to Figure 2d The spatial distribution pattern of key ecological zones in the Yellow River Basin from 2005 to 2020 is shown; Figure 3 Shows the classification of key ecological zones in the Yellow River Basin.

[0244] Furthermore, the above identification results were superimposed and compared with the most important nature reserves ECAs provided by the Ministry of Ecology and Environment. Specifically, the scope of ECAs includes priority areas for biodiversity conservation PABCs and national nature reserves NNRs.

[0245] The results of the above embodiment show that 23 of the 28 existing ECAs are within the range of CEAs identified in the above embodiment, such as Figures 4a to 4d As shown in , the recognition accuracy reaches 82%.

[0246] Except for the OHYZ in the West Ordos-Helan Mountain-Yinshan District, most of the most important nature reserves in other areas were identified as CEAs by the method in the above embodiment. The main reason for the low identification accuracy of the OHYZ in the West Ordos-Helan Mountain-Yinshan District is that the OHYZ mainly focuses on protecting desert ecosystems and geological relics, such as the Etuoke Dinosaur Relics and Fossils National Nature Reserve and the West Ordos National Nature Reserve in Inner Mongolia. Since the services and functions of desert ecosystems are not considered in the above embodiment, some of the most important nature reserves in the OHYZ are not identified as CEAs by the above embodiment.

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

[0248] The present invention has been described above in conjunction with preferred embodiments, but these embodiments are only exemplary and serve only as an illustration. On this basis, the present invention may be subjected to a variety of substitutions and improvements, all of which fall within the scope of protection of the present invention.

Claims

1. A method for identifying key ecological areas based on ecosystem multifunctionality and stability, characterized in that: The method comprises the following steps: Step 1: Obtain multi-source remote sensing datasets 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 within N years, and extract the spatial distribution pattern of ecosystem stability in the target area within 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.

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 2, the ecosystem multifunctionality index of the target area in the pixel scale for the spatial distribution pattern of ecosystem multifunctionality within N years is obtained by the following formulas (I) and (II): EMFI k = Mean(EMFI kj ) (1) Among them, EMFI k represents the ecosystem multifunctionality index in the kth year; EMFI kj represents the ecosystem multifunctionality index at the jth threshold level in the kth year; 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 ith 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 r i The value is 1, otherwise r i The value is 0; t ij represents the threshold of the ith ecosystem service under the threshold level; Mean(x) represents the mean function.

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 service supply stability index of the spatial pattern of ecosystem service supply stability in the target area within N years at the pixel scale is obtained by the following formula (III): Among them, ESSSI represents the ecosystem service supply stability index within N years; f(x) represents Jenks classification function; n represents the quantity of ecosystem service supply 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 in N years; SD pq It represents the variance of the qth ecological service supply on the pth grid within N years.

5. 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 spatial pattern of ecosystem resilience in the target area within N years at the pixel scale is obtained by the following formula (IV): ERI=f(Kendall τ(AR(1,NDVI),t))(4) Among them, ERI represents the ecosystem resilience index within N hands; f(x) represents 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 autoregression coefficient of the N-year NDVI dataset with a lag time of 1 month.

6. 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 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 (V): 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 Indicates the weight of ERI.

7. The method for identifying key ecological areas based on ecosystem multifunctionality and stability according to claim 1, characterized in that: In step 4, the spatial distribution pattern of the key ecological zones of the ecosystem in the target area within N years is obtained by the following equations (VI) and (VII) at the pixel scale: CEAs k = Top 30% (CEAI k )(VII) 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 the kth year; ESI represents the ecosystem stability index within N years; represents the normalization function; ω EMFI Represents EMFI k The weight of ω ESI Indicates the weight of ESI.

8. 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 areas, the key ecological area level is obtained at the pixel scale by 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.

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

10. The method for identifying key ecological areas based on ecosystem multifunctionality and stability according to claim 3, 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