Ecological protection area identification method based on ecological system service supply and demand mismatching

Abstract

The invention discloses an ecological protection area identification method based on ecological system service supply and demand mismatching, which comprises the following steps: measuring and drawing each ecological system service of a supply side and a demand side to obtain a supply horizontal spatial distribution map and a demand horizontal spatial distribution map; making a supply-demand aggregation graph by adopting GeoDa software, calculating by adopting ArcGIS software to obtain a comprehensive supply-demand ratio index, and constructing a comprehensive supply-demand ratio index spatial distribution graph; three main driving factors of the research area are screened out based on an OPGD model, and a driving space heterogeneity distribution diagram of the three driving factors is obtained by adopting an MGWR model; based on each ecological system service supply and demand aggregation graph, performing equal-weight integration by using a weighted stacking tool of ArcGIS software to generate a summary grid; then, based on the comprehensive supply-demand ratio index spatial distribution map, dividing the comprehensive supply-demand ratio index spatial distribution map into three clear areas by adopting a natural discontinuity point grading method; and finally, carrying out spatial integration by adopting a K-mode clustering analysis method in combination with the heterogeneity distribution diagram of the driving space, and completing the recognition of the ecological protection area.

Description

Technical Field

[0001] This invention belongs to the field of ecological space planning technology, specifically relating to a method for identifying ecological protection zones based on the mismatch between the supply and demand of ecosystem services. Background Technology

[0002] With the deepening of ecological civilization construction and the increasing demand for refined governance of national land space, the study of the supply and demand balance of ecosystem services, as a key link in maintaining regional ecological security and human well-being, has become a cutting-edge hot topic in geography, ecology, and resource and environmental management. However, the analysis and regulation of supply and demand mismatch of ecosystem services, especially in the study of complex systems across administrative regions, still faces significant technical bottlenecks.

[0003] Currently, existing technologies and research largely focus on single or a few key ecosystem services (such as the "water-carbon-food" paradigm). Static assessments of ecosystem service supply or demand at a single static scale (such as county or province) primarily rely on traditional statistical models (such as multiple regression and cluster analysis) and basic spatial analysis tools. These methods have fundamental limitations in analyzing the dynamics of ecosystem service supply and demand, cascading effects, and driving mechanisms. Specifically:

[0004] 1. Neglecting the cascade effect of "water-carbon-soil": Current technologies rarely characterize the cascade effect formed by ecosystem services through the "water-carbon-soil" functional chain (such as the regulatory role of soil structure on water conservation and carbon sequestration efficiency). The traditional "water-carbon-food" paradigm often simplifies soil conservation (SR) services as a passive support carrier for food production, neglecting the key "valve" role of soil conservation (SR) as an active regulatory factor in water supply (WY) and carbon sequestration (CS) and its cross-scale feedback mechanism.

[0005] 2. Inadequate Analysis of Driving Mechanisms: Traditional methods (such as ordinary least squares regression and simple geographic detectors) suffer from biases when analyzing the complex nonlinear and spatially heterogeneous impacts of multidimensional "social-ecological" driving factors (such as topography, land use change, urbanization intensity, and policies) on supply-demand mismatch. They struggle to handle multi-factor synergistic effects, spatiotemporal lag effects, and significant differences in driving forces at different spatial scales (such as county-level areas, small watersheds, and pan-plots), leading to inaccurate attribution and affecting the targeting of ecological protection zone strategies.

[0006] 3. The identification criteria for ecological protection zones are singular and their spatial adaptability is weak: Existing identification and zoning management of protected areas based on the mismatch between ecosystem service supply and demand often rely on a single service or a singular demand-side consideration. They fail to effectively integrate the multi-dimensional dynamic coupling relationship of "spatial correlation characteristics (synergy / trade-off) - supply and demand balance (surplus / deficit) - dominant driving mechanism," and also lack a systematic consideration of cross-regional scale effects at multiple levels (county-level administrative units and grid-based natural units). This results in blurred spatial boundaries for identified ecological protection zones (e.g., core conflict areas are masked at the county level) and insufficient adaptability, making it difficult to effectively support the need for refined collaborative governance based on natural geographical units and breaking through administrative barriers. Summary of the Invention

[0007] To address the aforementioned shortcomings of existing technologies, the present invention aims to provide a method for identifying ecological protection zones based on the mismatch between the supply and demand of ecosystem services. This invention can accurately identify ecological protection zones, thereby providing a scientific basis and technical support for the collaborative governance of ecological regions across administrative regions.

[0008] The technical solution of this invention is implemented as follows:

[0009] A method for identifying ecological protected areas based on the mismatch between the supply and demand of ecosystem services includes the following steps:

[0010] Step S1: Identify the key ecosystem service types in the study area;

[0011] Step S2: Based on the key ecosystem service types, determine the evaluation method for each ecosystem service type, find relevant data according to the corresponding evaluation method, measure and map the various ecosystem services on the supply and demand sides, and obtain spatial distribution maps of the supply level and demand level of each ecosystem service.

[0012] Step S3: Based on the spatial distribution maps of the supply and demand levels of various ecosystem services obtained in Step S2, GeoDa software is used to create supply and demand cluster maps of various ecosystem services to obtain the spatial relationship between the supply and demand of various ecosystem services. At the same time, ArcGIS software is used to calculate the comprehensive supply and demand ratio index and construct a spatial distribution map of the comprehensive supply and demand ratio index to measure the overall spatiotemporal surplus-deficit evolution pattern of supply and demand in the study area.

[0013] Step S4: Based on the OPGD model, three main driving factors in the study area were selected, and then the MGWR model was used to obtain the spatial heterogeneity distribution map of the three driving factors.

[0014] Step S5: Based on the ecosystem service supply and demand aggregation maps obtained in Step S3, use the weighted overlay tool in ArcGIS software to perform equal-weighted aggregation, generating a summary raster that reflects the spatial coordination / trade-off characteristics of ecosystem service supply and demand; then, based on the spatial distribution map of the comprehensive supply-demand ratio index obtained in Step S3, use the natural discontinuity grading method to divide it into three distinct regions; finally, combine the driving spatial heterogeneity distribution map obtained in Step S4, and use the K-modes clustering analysis method to perform spatial integration, completing the identification of ecological protection zones.

[0015] Furthermore, the key ecosystem services are water supply, carbon sequestration, and soil conservation.

[0016] Furthermore, in step S2, the three ecosystem services of water supply, carbon sequestration and soil conservation on the supply side are measured and mapped based on the INVEST model, CASA model and RUSLE model, respectively, to obtain spatial distribution maps of the supply levels of each ecosystem service.

[0017] Furthermore, the supply-side evaluation formula for water supply is as follows:

[0018]

[0019]

[0020] In the formula, WY x The water supply measurement value for the grid cell x supply side, in mm; AET x P represents the actual annual evapotranspiration of grid cell x, in mm. x R represents the annual precipitation of grid cell x, in mm; x Z is the drying index of grid cell x, dimensionless; Z is the Zhang coefficient, dimensionless; w x AWCx is the ratio of annual water requirement to precipitation for vegetation; AWCx is the available water for plants on grid cell x; ET ox Potential evapotranspiration in grid cell x, mm; K x This represents the vegetation coefficient.

[0021] The formula for evaluating the supply side of carbon sequestration is as follows:

[0022] NPP x =APAR x ×ε x

[0023] APAR x =TSOL·FPAR·r

[0024]

[0025] In the formula, NPPx This represents the carbon retention measure on the supply side of grid cell x, in gC / m. 3 APAR x and ε x These represent the photosynthetically active radiation and actual light energy utilization rate at grid cell x, respectively. TSOL represents the total solar radiation irradiated to the Earth's surface, in MJ·m³. -2 FPAR represents the amount of solar radiation absorbed by the vegetation layer for photosynthesis, measured in MJ·m³. -2 ; r represents the ratio of solar radiation absorbed by vegetation for photosynthesis to total solar radiation, dimensionless; SI represents the ratio vegetation index, ranging from 4.14 to 6.17; NDVI is the normalized difference vegetation index.

[0026] The formula for evaluating the supply side of soil conservation is as follows:

[0027] A c =A p -A r =R×K×L×S×(1-CP)

[0028] In the formula, A c For the soil retention measurement value on the supply side, t / hm 2 A p Potential erosion, t / hm 2 A r This represents the actual erosion rate, in t / hm². 2 R represents the rainfall erosion factor, in MJ·mm·ha. -1 ·h -1 ·a -1 K is the erosion factor t·ha·MJ -1 ·mm -1 L is the slope length factor; S is the slope gradient factor; C is the vegetation cover factor; P is the water conservation factor. L, S, C, and P are all dimensionless factors.

[0029] Furthermore, based on ArcGIS software, the two ecosystem services of water supply and carbon sequestration on the demand side were measured and mapped to obtain spatial distribution maps of the demand levels for water supply and carbon sequestration; based on the RUSLE model, the ecosystem service of soil conservation on the demand side was measured and mapped to obtain spatial distribution maps of the demand levels for soil conservation.

[0030] Furthermore, the formula for evaluating the supply and demand sides of water resources is as follows:

[0031] D sum =D a ×p pop +D b ×GDP+D c ×AGR

[0032] In the formula, D sum Water supply measurement value on the demand side, mm; D a per capita water demand (m) 3 ;p pop For population density data, people / hm 2 ;D b Water consumption per 10,000 yuan of GDP, in meters; GDP is GDP data, 10,000 yuan / hm². 2 ;D c The amount of irrigation water per unit area of ​​cultivated land, in m 3 AGR represents farmland raster data.

[0033] The formula for evaluating the demand side of carbon sequestration is as follows:

[0034] D cs,x =C e.x ×p x ×C t

[0035] In the formula, D cs,x The carbon retention measure at grid cell x on the demand side is gC / m. 3 C e,x p represents the per capita energy consumption at grid cell x, in 10,000 standard coal equivalents. x C represents the population density at grid cell x; x The carbon emission conversion rate of energy consumption is 0.6739.

[0036] The soil conservation demand-side evaluation formula is as follows:

[0037] A p =R×K×L×S×P

[0038] Furthermore, based on the measures of various ecosystem services obtained in steps S1 and S2, a comprehensive supply-demand ratio index is calculated to measure the overall spatiotemporal surplus-deficit evolution of supply and demand in the study area. The comprehensive supply-demand ratio index is calculated using the following formula:

[0039]

[0040] In the formula: CESDR is the composite supply-demand ratio index. A positive CESDR indicates oversupply, a CESDR of 0 indicates supply-demand balance, and a negative CESDR indicates undersupply; S i D is a measure of the supply side of the i-th ecosystem service; i S is a measure of the demand side of the i-th ecosystem service; imax D represents the maximum measure of the supply side of the i-th ecosystem service; imaxis the maximum measure of the demand side of the i-th ecosystem service; n is the total number of ecosystem services.

[0041] Furthermore, the driving factors include land use type, soil type, elevation, temperature, solar radiation, precipitation, potential evapotranspiration, vegetation index, vegetation available water content, GDP, and population density.

[0042] Furthermore, step S6 is included: after identifying the ecological protection zone, if there is a conflict with the local economic development plan, the weights of various ecosystem services are adjusted, and then the weighted overlay tool of ArcGIS software is used to synthesize the data to generate a summary raster that reflects the spatial coordination / trade-off characteristics of ecosystem service supply and demand. Then, step S5 is repeated to complete the identification of the ecological protection zone.

[0043] Compared with the prior art, the present invention has the following beneficial effects:

[0044] 1. This invention takes water supply, carbon sequestration, and soil conservation as key ecosystem services, constructs a "water-carbon-soil" system model, quantifies the cascade regulatory effect of soil conservation on water supply and carbon sequestration, breaks through the limitation of "emphasizing supply while neglecting flow", realizes a comprehensive analysis of the supply and demand relationship of ecosystem services, can significantly improve the spatial continuity and comprehensiveness of ecological protection areas, thereby effectively solving the spatial fragmentation problem caused by traditional methods and enhancing the overall protection effectiveness of ecosystems.

[0045] 2. This invention integrates the OPGD model and the MGWR model, and incorporates the spatiotemporal heterogeneous contributions of multidimensional "social-ecological" factors, significantly improving the accuracy of driving mechanism analysis. This enables the generation of a more precise library of differentiated regulation strategies, achieving targeted protection adaptation for ecological protection zones and effectively addressing complex and ever-changing ecological and environmental challenges.

[0046] 3. This invention proposes a "ternary coupling weighting method," which integrates the multi-dimensional dynamic coupling relationship of "spatial correlation characteristics (synergy / trade-off) - supply and demand balance state (surplus / deficit) - dominant driving mechanism" to achieve the scientific delineation of imbalance hotspots and realize the dynamic identification of ecological protection zones based on multi-level (county-level administrative units and grid natural units) cross-regional scale effects. Furthermore, the method described in this invention not only has strong coupling but can also be effectively transferred to other urban agglomeration scales to achieve effective ecological protection zone delineation.

[0047] 4. This invention includes a comprehensive zoning method encompassing three dimensions: spatial correlation characteristics, supply and demand balance, and dominant driving factors. Combined with a two-level management and grading strategy of "county-level leadership + grid correction," it can provide strong support for the coordinated advancement of ecological protection and economic development. Attached Figure Description

[0048] Figure 1- Technical roadmap of the present invention.

[0049] Figure 2 - Spatial distribution map of water supply levels.

[0050] Figure 3 - Spatial distribution of carbon sequestration supply levels.

[0051] Figure 4 - Spatial distribution map of soil retention supply levels.

[0052] Figure 5 - Spatial distribution map of water supply demand levels.

[0053] Figure 6 - Spatial distribution map of carbon sequestration demand levels.

[0054] Figure 7 - Spatial distribution map of soil conservation needs.

[0055] Figure 8 - Supply and demand aggregation diagram of various ecosystem services.

[0056] Figure 9 - Spatial distribution map of the comprehensive supply-demand ratio index.

[0057] Figure 10 -Results of single-factor detection of the comprehensive supply-demand ratio index.

[0058] Figure 11 - Results of dual-factor interactive detection of the comprehensive supply-demand ratio index.

[0059] Figure 12 - Distribution of heterogeneity in the driving space of driving factors.

[0060] Figure 13 - Mapping of ecological protection zones at the county level.

[0061] Figure 14 - Grid-scale ecological protection zone delineation map. Detailed Implementation

[0062] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments.

[0063] This embodiment takes the Chengdu-Chongqing Twin-City Economic Circle as the study area, and then identifies ecological protection zones based on the mismatch between the supply and demand of ecosystem services. The main data sources are shown in Table 1, where the data is from the study area in 2020. Figures 2-9 The year 2020 in the text refers to the year.

[0064] Table 1. Main Data Sources

[0065]

[0066] To ensure the accuracy and reliability of the data, the following special notes are provided:

[0067] (1) GDP is calibrated at the county level based on the China Urban Statistical Yearbook. POP is calibrated at the county level based on population census data.

[0068] (2) Monthly NDVI data and solar radiation data were extracted through the GEE platform and verified according to Xu Xinliang. Data set of NDVI and EVI250m every 16 days in China. Resource and Environmental Science Data Registration and Publication System (http: / / www.resdc.cn / DOI), 2024.DOI:10.12078 / 2024072902 and China Meteorological Data Network (https: / / data.cma.cn).

[0069] The method for identifying ecological protected areas based on the mismatch between the supply and demand of ecosystem services includes the following steps:

[0070] (1) Based on the inherent differences among the three ecosystem services of water supply (WY), carbon sequestration (CS) and soil conservation (SR), the INVEST model, CASA model and RUSLE model were used to measure and map the three ecosystem services of water supply, carbon sequestration and soil conservation on the supply side, respectively, and spatial distribution maps of the supply level of each ecosystem service were obtained.

[0071] The formula for evaluating water supply on the supply side is as follows:

[0072]

[0073]

[0074] In the formula, WY x The water supply quantity for grid cell x, in mm, is the measured value of the water supply on the supply side of grid cell x; AET x P represents the actual annual evapotranspiration of grid cell x, in mm. x R represents the annual precipitation of grid cell x, in mm; x Z is the drying index of grid cell x, dimensionless; Z is the Zhang coefficient, dimensionless; w x AWCx is the ratio of annual water requirement to precipitation for vegetation; AWCx is the available water for plants on grid cell x; ET ox Potential evapotranspiration in grid cell x, mm; K x The vegetation coefficient represents the ratio of plant evapotranspiration to potential evapotranspiration at different developmental stages. The resulting spatial distribution map of water supply levels is shown below. Figure 2As shown in the figure, the water supply service volume ranges from [0, 1513.3]. The western and northwestern Minshan-Qionglai-Liangshan region (Ya'an City, southwestern Leshan City, and northwestern Chengdu, Deyang, and Mianyang) and the strip-shaped mountain ranges of the Huaying Mountains, Tongluo Mountains, and Mingyue Mountains in northeastern Chongqing are characterized by low ES values. The Chengdu urban area and Chongqing main urban area, dominated by construction land, are characterized by high values. Furthermore, the eastern edge of the Minshan-Qionglai-Liangshan region also belongs to a high-value area, with the maximum water supply depth reaching 1513.3 mm.

[0075] The formula for evaluating the supply side of carbon sequestration is as follows:

[0076] NPP x =APAR x ×ε x

[0077] APAR x =TSOL·FPAR·r

[0078]

[0079] In the formula, NPP x Represents the net primary productivity at grid cell x, in gC / m³. 3 This refers to the carbon retention measure on the supply side of the grid cell x; APAR x and ε x These represent the photosynthetically active radiation and actual light energy utilization rate at grid cell x, respectively. TSOL represents the total solar radiation irradiated to the Earth's surface, in MJ·m³. -2 FPAR represents the amount of solar radiation absorbed by the vegetation layer for photosynthesis, measured in MJ·m³. -2 ; r represents the ratio of solar radiation absorbed by vegetation for photosynthesis to total solar radiation, dimensionless; SI represents the ratio vegetation index, ranging from 4.14 to 6.17. Inputting the vegetation index NDVI into the CACA model will yield the maximum and minimum values ​​of SI; NDVI is the normalized difference vegetation index. The resulting spatial distribution map of carbon sequestration supply levels is shown below. Figure 3 As shown in the figure, the carbon sequestration service supply ranges from [0, 217.52], and spatially it generally exhibits a distribution pattern of "high in the periphery and low in the two cores". The strip-shaped mountain ranges consisting of the Minshan-Qionglai-Liangshan region (Ya'an City, southwestern Leshan City, and northwestern Chengdu, Deyang City, and Mianyang City) in the west and northwest, as well as the Huaying Mountains, Tongluo Mountains, and Mingyue Mountains in northeastern Chongqing, are characterized by high-value ES distribution. The urban area of ​​Chengdu and the main urban area of ​​Chongqing, which are mainly composed of construction land, are characterized by low-value areas, and the area of ​​both regions has expanded over time.

[0080] The formula for evaluating the supply side of soil conservation is as follows:

[0081] A c =A p -A r =R×K×L×S×(1-CP)

[0082] In the formula, A c Soil retention capacity, t / hm 2 That is, the soil retention measurement value on the supply side; A p Potential erosion, t / hm 2 A r This represents the actual erosion rate, in t / hm². 2 R represents the rainfall erosion factor, in MJ·mm·ha. -1 ·h -1 ·a -1 K is the erosion factor t·ha·MJ -1 ·mm -1 L represents slope length; S represents slope gradient; C represents vegetation cover; and P represents water conservation. L, S, C, and P are all dimensionless factors. The resulting spatial distribution map of soil conservation supply levels is shown below. Figure 4 As shown in the figure, the soil conservation service supply ranges from [0, 1613.38], exhibiting a spatial distribution pattern of "high in the periphery and low in the interior." The strip-shaped mountain ranges forming the western and northwestern Minshan-Qionglai-Liangshan region (Ya'an City, southwestern Leshan City, and northwestern Chengdu, Deyang City, and Mianyang City) and the northeastern Chongqing region, consisting of the Huaying Mountains, Tongluo Mountains, and Mingyue Mountains, are characterized by high ES values.

[0083] (2) Based on ArcGIS software, the two ecosystem services of water supply and carbon sequestration on the demand side were measured and mapped to obtain spatial distribution maps of the demand levels of water supply and carbon sequestration; based on the RUSLE model, the ecosystem service of soil conservation on the demand side was measured and mapped to obtain spatial distribution maps of the demand levels of soil conservation.

[0084] The water supply and demand assessment considers the actual consumption of freshwater resources, that is, people's demand for water resources, including agricultural production water demand, industrial water demand, and domestic water demand, as shown in the following formula:

[0085] D sum =D a ×p pop +D b ×GDP+D c ×AGR

[0086] In the formula, D sum The demand for water supply services, expressed in mm, is the demand-side measure of water supply; D a per capita water demand (m)3 ;p pop For population density data, people / hm 2 ;D b Water consumption per 10,000 yuan of GDP, in meters; GDP is GDP data, 10,000 yuan / hm². 2 ;D c The amount of irrigation water per unit area of ​​cultivated land, in m 3 AGR represents cultivated land raster data. The resulting spatial distribution map of water supply and demand is shown below. Figure 5 As shown in the figure, the demand for water supply services ranges from [0, 5611], exhibiting a spatial distribution pattern of "high demand in the two core cities, moderate demand in the central region, and low demand in the periphery." The dual-core cities of Chengdu and Chongqing have a significantly higher demand for water supply services (ES) than other prefecture-level cities or counties. This is because these areas are densely populated and have a higher water demand, with a maximum water depth of 5611 mm. In contrast, sparsely populated areas, such as Ya'an in the west and Luzhou in the south, have relatively lower demand. Notably, the areas with lower demand are characterized by forest land having lower demand than cultivated land, because forest land has fewer crops and therefore lower water requirements.

[0087] The demand-side assessment of carbon sequestration considers regional carbon emissions as a conservative estimate of the demand for carbon sequestration services, as shown in the following formula:

[0088] D cs,x =C e.x ×p x ×C t

[0089] In the formula, D cs,x The carbon retention measure at grid cell x on the demand side is gC / m. 3 C e,x p represents the per capita energy consumption at grid cell x, in 10,000 standard coal equivalents. x C represents the population density at grid cell x; x The carbon emission conversion rate of energy consumption is taken as 0.6739. The spatial distribution map of carbon sequestration demand is shown below. Figure 6 As shown in the figure, the demand for carbon sequestration services ranges from [0, 632.557], exhibiting a spatial distribution pattern of "high demand in the two core cities, moderate demand in the central region, and low demand in the periphery." The demand for ES in the dual-core cities of Chengdu and Chongqing is significantly higher than in other prefecture-level cities or counties. This is because these areas are densely populated, with carbon demand reaching 1182.81 gC / m³. 2 In sparsely populated areas, such as Ya'an in the west and Luzhou in the south, demand is relatively low.

[0090] Soil conservation demand-side assessment considers actual soil erosion as the amount of ecosystem services that humans expect to manage and obtain, characterizing the demand for soil conservation services by the actual amount of soil erosion. The formula is as follows:

[0091] A p =R×K×L×S×P

[0092] In the formula, A p Potential erosion, t / hm 2 , which is the soil retention measure on the demand side; R is the rainfall erosion factor, MJ·mm·ha -1 ·h -1 ·a -1 K is the erosion factor, t·ha·MJ -1 ·mm -1 L represents the slope length factor; S represents the slope gradient factor; and P represents the water conservation factor. L, S, and P are all dimensionless factors. The spatial distribution map of soil conservation requirements is shown below. Figure 7 As shown in the figure, the demand for soil conservation services (SR) is in the range of [0, 1702.77], exhibiting a similar "high in the periphery, low in the interior" distribution pattern to the supply side. Specifically, the high-value areas for SR are the strip-shaped mountain ranges formed by the Minshan-Qionglai-Liangshan region (Ya'an City, southwestern Leshan City, and northwestern Chengdu, Deyang City, and Mianyang City) in the west and northwest, as well as the Huaying Mountains, Tongluo Mountains, and Mingyue Mountains in northeastern Chongqing.

[0093] (3) Using GeoDa software, an LISA (Local Indicators of Spatial Association) clustering map was created to demonstrate the degree of clustering between the supply-side and demand-side measures of the three ecosystem services. The resulting ecosystem service supply-demand clustering map is shown below. Figure 8 As shown in the figure, the clusters are divided into High-High Cluster (HH), Low-Low Cluster (LL), High-Low Outlier (HL), Low-High Outlier (LH), and No Significant (N). This visually reflects the spatial relationship between the supply and demand of various ecosystem services. HH and LL clusters indicate a synergistic relationship, while HL and LH clusters indicate a trade-off relationship.

[0094] (4) Based on the measurement values ​​of various ecosystem services obtained in steps (1) and (2), the supply-demand ratios of water supply, carbon sequestration, and soil conservation are calculated to characterize the balance between the supply and demand of each ecosystem service. A comprehensive supply-demand ratio index is constructed to reflect the spatiotemporal surplus-deficit evolution of the overall supply and demand in the Hengliang study area. The formula for calculating the comprehensive supply-demand ratio index is as follows:

[0095]

[0096] In the formula: CESDR is the composite supply-demand ratio index. A positive CESDR indicates oversupply, a CESDR of 0 indicates supply-demand balance, and a negative CESDR indicates undersupply; S i Let D be the supply of the i-th ecosystem service, i.e., the measure of the supply side of the i-th ecosystem service; i Let S be the demand for the i-th ecosystem service, which is a measure of the supply side of the i-th ecosystem service; imax D represents the maximum supply of the i-th ecosystem service; imax Let n be the maximum demand for the i-th ecosystem service; n is the total number of ecosystem services. In this embodiment, n = 3, and i is 1, 2, or 3.

[0097] The spatial distribution diagram of the comprehensive supply-demand ratio index obtained in this embodiment is as follows: Figure 9 As shown in the figure, CESDR accounts for 0.5% in [-0.85, -0.2); 0.6% in [-0.2, -0.1); 0.9% in [-0.1, -0.05); 1.3% in [-0.05, 0); 7% in [0, 0.05); 34.9% in [0.05, 0.1); 46.3% in [0.1, 0.2); and 8.6% in [0.2, 0.65].

[0098] (5) Three main driving factors of the study area were selected based on the OPGD model.

[0099] Traditional methods require manual classification of variable categories and quantities, resulting in strong subjectivity and poor computational accuracy. OPGD, however, can optimize the discretization of continuous variables through various statistical methods (equal intervals, natural intervals, quantile intervals, geometric intervals, and standard deviation intervals) to obtain the optimal q-value, thus more accurately reflecting driving forces. Therefore, OPGD is used to analyze the influencing factors of CESDR. The factor detector and interaction detector of OPGD are used to analyze the influencing factors of CESDR. The factor detector detects the degree to which the independent variable X (influencing factor) explains the spatial differentiation of Y; the interaction detector is used to identify the interaction between different risk factors Xs, that is, to assess whether the explanatory power of factors X1 and X2 on the dependent variable Y is enhanced, weakened, linear, or nonlinear. The formula for calculating the q-value is:

[0100]

[0101] In the formula, q measures the explanatory power of the influencing factors on Y, with a range of [0,1]. A larger q value indicates stronger explanatory power from each influencing factor, meaning a more significant impact on the spatial distribution of CESDR. L represents the stratification of Y or X; N... h σ represents the number of elements in layer h and the entire region; h 2SST represents the variance of the Y values ​​for layer h and the entire region; SST represents the total variance of the entire region.

[0102] Selection of Ecosystem (ES) Drivers. Topography is an important factor influencing the spatial differentiation of ES. It can affect the distribution of water and heat energy, thereby influencing the distribution patterns of water production and soil erosion. Climate conditions affect vegetation distribution, providing a good foundation for carbon sequestration. Population size and GDP, which indirectly reflect the intensity of human activities, also interfere with the supply and demand differences of ES. Based on the above factors, 11 driving factors were selected from two dimensions: natural environment and socio-economic (Table 2). For continuous data, optimal classification and determination of the number of classifications are required. The optimal parameter geographic exploration analysis is carried out using the comprehensive supply-demand ratio (CESDR) calculated in step (4), aiming to explain the main driving forces of the overall ecological function of the study area.

[0103] Table 2. Explanation of Driving Factors and Grille Data

[0104]

[0105] At the optimal research scale, existing studies typically compare the 90th quantile of the q-values ​​of all driving factors at each spatial scale. The optimal spatial scale is determined when the quantile reaches its maximum. For urban agglomeration scales, spatial scale units are generally set at 3, 5, 8, 10, and 15 km to determine the optimal scale. For discretization methods, equal intervals, natural discontinuities, quantiles, geometric intervals, and standard deviations are selected, with the number of classifications set to 5-10. The parameter with the largest q-value is selected as the optimal parameter for the geographic detector analysis. Taking the Chengdu-Chongqing Twin-City Economic Circle as an example, the 90th quantile reaches its peak at 10 km according to the detection results. Therefore, among the five grid types, the 10 km grid better reflects the impact of driving factors on CESDR. Secondly, regarding the optimal discretization method and the optimal number of partitions, for GDP (X3) and POP (X4), q is maximized when the classification method is geometric interval classification and the number of partitions is 10. Therefore, the geometric interval classification method is used to divide GDP (X3) and POP (X4) into 10 categories. Similarly, the standard deviation method is used to divide solar radiation (X7) and precipitation (X8) into 9 categories; the natural discontinuity method is used to divide air temperature (X6) and NDVI (X10) into 10 and 8 categories, respectively; and the quantile method is used to divide DEM (X5), potential evapotranspiration (X9), and PAWC (X11) into 9, 9, and 10 categories, respectively. Based on the determination of the above optimal parameters, the results of single-factor and two-factor interactive CESDR detection are further obtained. The single-factor interactive results and multi-factor interactive results are as follows: Figure 10 and Figure 11 As shown, where Figure 10The asterisk (*) indicates a significance level of 0.01. As shown in the figure, the single-factor ranking is: 0.57 (POP) > 0.55 (GDP) > 0.37 (NDVI) > 0.32 (Elevation) > 0.29 (Land Use Type) > 0.29 (Temperature) > 0.24 (Soil Type) > 0.20 (Potential Evapotranspiration) > 0.19 (PAWC) > 0.16 (Precipitation) > 0.15 (Solar Radiation). The main influencing factors are GDP (X3) and POP (X4), with driving forces of 0.55 and 0.57 respectively. NDVI (X10) is the second most influential factor, with a driving force of 0.37. This indicates that human factors are a key factor affecting the sustainable development of ES. This is further supported by the results of the two-factor interaction detection. The combined effect of different factors enhanced the explanatory power of CESDR, indicating that the influence of ES is often the result of the synergistic effect of multiple factors. The synergistic driving force of GDP(X3)∩NDVI(X10) (q=0.67) and POP(X4)∩NDVI(X10) (q=0.70) was the most prominent, which means that promoting the coordination of economic development and ecological protection is crucial to regional ecological security. Therefore, these three main driving factors were selected as the primary driving factors.

[0106] (6) The heterogeneity distribution of the driving space of the three driving factors was obtained by using the MGWR model.

[0107] After identifying the main driving factors, it is also necessary to understand their spatial effects. Traditional geographically weighted regression (GWR) models, due to the constant bandwidth settings for different explanatory variables, largely limit the variation of spatially non-stationary relationships. To address the limitations of GWR, the MGWR model is used here, which effectively solves its shortcomings. The calculation formula for the MGWR model is as follows:

[0108]

[0109] In the formula, (u i ,v i ) represents the position of the i-th grid cell; Y i and X ij These are the distribution density Y and the influence factor X, respectively. j in (u i v i The observation at ); β bwj β is the regression coefficient of the j-th influencing factor under the optimal bandwidth. bw0 For (u i ,v i The intercept at ); b wj and b w0 The optimal bandwidths used for the j-th influence factor and intercept are respectively; E i Let be the independently distributed random residual of the i-th grid cell. The resulting heterogeneity distribution map of the driving space is as follows: Figure 12 As shown.

[0110] (7) Based on the ecosystem service supply and demand aggregation maps obtained in step (3), use the weighted overlay tool of ArcGIS software to perform equal-weighted aggregation. During the overlay process, for areas with the same proportion of various ecosystem service aggregation types, set priority principles according to the degree of ecological security threat (LH>HH=LL>HL>N), and retain the original spatial correlation zoning structure (coordination zone: HH / LL; trade-off zone: HL / LH) to generate a summary raster reflecting the overall spatial coordination / trade-off characteristics of ecosystem service supply and demand. Then, based on the spatial distribution map of the comprehensive supply and demand ratio index obtained in step (4), use the natural discontinuity grading method to divide it into three distinct areas: the deficit core area with CESDR<-0.1 (representing severe supply). The transition buffer zone with a CESDR of -0.1≤CESDR≤0.1 (representing a basic balance or slight fluctuation in supply and demand) and the surplus protection zone with a CESDR>0.1 (representing a relatively abundant supply) are defined. Finally, the spatial heterogeneity distribution map of the three driving factors, GDP (economic development), POP (population pressure) and NDVI (vegetation cover), obtained in step (5), is used to identify three types of dominant regions: the social pressure-dominated region (strong GDP and POP effects, weak NDVI effects), the multi-driving coupling region (moderate GDP, POP, and NDVI effects), and the ecological restoration sensitive region (strong NDVI effects, weak GDP and POP effects). Based on the analysis results from the three dimensions mentioned above (spatial correlation, supply and demand status, and driving mechanism), K-mode clustering was used for spatial integration. This was then validated and adjusted in conjunction with the local economic development plan of the study area (CCEC) (specifically by adjusting the weights of various ecosystem services). Ultimately, three types of ecological management functional zoning schemes were formed: core restoration areas, resilient regulation areas, and ecological priority conservation areas (ecological protection areas), thus completing the identification of ecological protection areas. Based on the ecological management functional zoning results, a two-tiered management and grading strategy of "county-level leadership + grid correction" was generated, encompassing natural restoration, economic compensation, and institutional synergy. The resulting county-level and grid-level ecological protection area zoning maps are shown below. Figure 13 and Figure 14 As shown in the figure:

[0111] Core Restoration Area: This area comprises five districts / counties: Yuzhong, Jiangbei, Nan'an, Jiulongpo, and Dadukou. These are the core areas of Chongqing, densely populated and highly modernized, with land use primarily for construction, and are the economic centers of Chongqing. This region suffers from a severe supply-demand deficit in ecosystem services (ES), primarily characterized by a concentration of low-lying vegetation (LH) and is dominated by social pressures, indicating a serious disruption to the supply and demand mechanisms of ecosystem services.

[0112] The resilient regulation zone comprises 13 districts and counties in Chongqing, including Banan, Shapingba, and Yubei districts, and 8 cities in Sichuan, including Chengdu, Deyang, and Ziyang. These areas are predominantly low hills with arable land as the main land use type. They are significantly impacted by human activities, resulting in low habitat quality, susceptibility to precipitation erosion, and a high risk of soil erosion. The region exhibits relatively low CESDR values, primarily acting as a transitional buffer zone on the edge of equilibrium, and is mainly characterized by high-risk (HH) accumulation. This indicates significant supply and demand pressures on ecosystem services in this region.

[0113] Ecological Priority Conservation Area (Ecological Protection Area): This region encompasses 18 cities and counties surrounding the CCEC, including 7 cities in Sichuan Province (such as Mianyang and Ya'an) and 11 districts and counties in Chongqing (such as Jiangjin, Qijiang, and Kaizhou). It is influenced by natural climate, high altitude, good vegetation productivity, abundant soil organic carbon reserves, and a strong ecological foundation. These areas have a high supply and demand for ecosystem services (ES), all belonging to surplus protection areas, and exhibiting clusters of both high-yield (HL) and low-yield (LL) ecosystem services. This indicates an ample supply of ecosystem services.

[0114] Compared to the county-level scale, the core restoration area accounts for a higher proportion (14.1%) at the grid-level. Furthermore, except for the five main urban districts of Chongqing, major urban areas in central China, as well as the urban area of ​​Chengdu, have all changed from flexible control areas to core restoration areas. This is because the proportion of urban core areas within administrative regions is relatively small, and from an administrative region scale, the overall ecological level is good. The proportion of ecological priority conservation areas decreased from 62.3% at the county-level scale to 28.5% at the grid-level (a decrease of 33.8%). This is because the study area is mainly surrounded by ecological barrier areas, encompassing the ecological barriers of the Minshan-Qionglai-Liangshan region, the Dalou Mountains, the Wuling Mountains, and the Daba Mountains in the west and northwest. These areas serve as the ecological baseline for each city and county, coordinating the impacts of urbanization. Therefore, at the administrative region scale, they are ecological priority protection areas, while the grid-level scale more accurately identifies these ecological priority protection areas.

[0115] Finally, it should be noted that the above embodiments of the present invention are merely illustrative examples and not intended to limit the implementation of the invention. Those skilled in the art can make other variations and modifications based on the above description. It is impossible to exhaustively list all possible implementations here. All obvious variations or modifications derived from the technical solutions of this invention are still within the scope of protection of this invention.

Claims

1. A method for identifying ecological protected areas based on the mismatch between the supply and demand of ecosystem services, characterized in that, Specifically, the following steps are included: Step S1: Identify the key ecosystem service types in the study area; Step S2: Based on the key ecosystem service types, determine the evaluation method for each ecosystem service type, find relevant data according to the corresponding evaluation method, measure and map the various ecosystem services on the supply and demand sides, and obtain spatial distribution maps of the supply level and demand level of each ecosystem service. Step S3: Based on the spatial distribution maps of the supply and demand levels of various ecosystem services obtained in Step S2, GeoDa software is used to create supply and demand cluster maps of various ecosystem services to obtain the spatial relationship between the supply and demand of various ecosystem services. At the same time, ArcGIS software is used to calculate the comprehensive supply and demand ratio index and construct a spatial distribution map of the comprehensive supply and demand ratio index to measure the overall spatiotemporal surplus-deficit evolution pattern of supply and demand in the study area. Step S4: Based on the OPGD model, three main driving factors in the study area were selected, and then the MGWR model was used to obtain the spatial heterogeneity distribution map of the three driving factors. Step S5: Based on the ecosystem service supply and demand aggregation maps obtained in Step S3, use the weighted overlay tool of ArcGIS software to perform equal-weighted aggregation and generate a summary raster that reflects the spatial coordination / trade-off characteristics of ecosystem service supply and demand. Then, based on the spatial distribution map of the comprehensive supply-demand ratio index obtained in step S3, it is divided into three distinct regions using the natural discontinuity grading method; finally, combined with the spatial heterogeneity distribution map of the driving space obtained in step S4, the K-modes clustering analysis method is used for spatial integration to complete the identification of the ecological protection zone.

2. The method for identifying ecological protected areas based on ecosystem service supply-demand mismatch according to claim 1, characterized in that, The key ecosystem services are water supply, carbon sequestration, and soil conservation.

3. The method for identifying ecological protected areas based on ecosystem service supply-demand mismatch according to claim 2, characterized in that, In step S2, the three ecosystem services of water supply, carbon sequestration and soil conservation on the supply side are measured and mapped based on the INVEST model, CASA model and RUSLE model, respectively, to obtain the spatial distribution map of the supply level of each ecosystem service.

4. The method for identifying ecological protected areas based on the mismatch between supply and demand of ecosystem services according to claim 3, characterized in that, The formula for evaluating water supply on the supply side is as follows: In the formula, WY x The water supply measurement value for the grid cell x supply side, in mm; AET x P represents the actual annual evapotranspiration of grid cell x, in mm. x R represents the annual precipitation of grid cell x, in mm; x Z is the drying index of grid cell x, dimensionless; Z is the Zhang coefficient, dimensionless; w x AWCx is the ratio of annual water requirement to precipitation for vegetation; AWCx is the available water for plants on grid cell x; ET ox Potential evapotranspiration in grid cell x, mm; K x The vegetation coefficient; The formula for evaluating the supply side of carbon sequestration is as follows: NPP x =FIRE x ×ε x MONKEYS x =TSOL FPAR r In the formula, NPP x This represents the carbon retention measure on the supply side of grid cell x, in gC / m. 3 APAR x and ε x These represent the photosynthetically active radiation and actual light energy utilization rate at grid cell x, respectively. TSOL represents the total solar radiation irradiated to the Earth's surface, in MJ·m³. -2 FPAR represents the amount of solar radiation absorbed by the vegetation layer for photosynthesis, measured in MJ·m³. -2 ; r represents the ratio of solar radiation absorbed by vegetation for photosynthesis to total solar radiation, dimensionless; SI represents the ratio vegetation index, ranging from 4.14 to 6.17; NDVI is the normalized difference vegetation index. The formula for evaluating the supply side of soil conservation is as follows: A c =A p -A r =R×K×L×S×(1-C-P) In the formula, A c For the soil retention measurement value on the supply side, t / hm 2 A p Potential erosion, t / hm 2 A r This represents the actual erosion rate, in t / hm². 2 R represents the rainfall erosion factor, in MJ·mm·ha. -1 ·h -1 ·a -1 K is the erosion factor t·ha·MJ -1 ·mm -1 L is the slope length factor; S is the slope gradient factor; C is the vegetation cover factor; P is the water conservation factor. L, S, C, and P are all dimensionless factors.

5. A method for identifying ecological protected areas based on ecosystem service supply-demand mismatch according to claim 2 or 3, characterized in that, Based on ArcGIS software, we measured and mapped the two ecosystem services on the demand side, water supply and carbon sequestration, to obtain spatial distribution maps of the demand levels for water supply and carbon sequestration. Based on the RUSLE model, the demand-side soil conservation ecosystem service is measured and mapped to obtain a spatial distribution map of the demand level for soil conservation.

6. The method for identifying ecological protected areas based on the mismatch between supply and demand of ecosystem services according to claim 5, characterized in that, The formula for evaluating the supply and demand of water resources is as follows: D sum =D a ×p pop +D b ×GDP+D c ×AGR In the formula, D sum Water supply measurement value on the demand side, mm; D a per capita water demand (m) 3 ;p pop For population density data, people / hm 2 ; D b Water consumption per 10,000 yuan of GDP, in meters; GDP is GDP data, 10,000 yuan / hm². 2 ; D c The amount of irrigation water per unit area of ​​cultivated land, in m 3 AGR represents farmland raster data. The formula for evaluating the demand side of carbon sequestration is as follows: D cs,x =C e.x ×p x ×C t In the formula, D cs,x The carbon retention measure at grid cell x on the demand side is gC / m. 3 C e,x p represents the per capita energy consumption at grid cell x, in 10,000 standard coal equivalents. x C represents the population density at grid cell x; x The carbon emission conversion rate of energy consumption is 0.6739. The soil conservation demand-side evaluation formula is as follows: A p =R×K×L×S×P。 7. The method for identifying ecological protected areas based on ecosystem service supply-demand mismatch according to claim 5, characterized in that, Based on the measures of various ecosystem services obtained in steps S1 and S2, the comprehensive supply-demand ratio index is calculated to measure the overall spatiotemporal surplus-deficit evolution of supply and demand in the study area. The comprehensive supply-demand ratio index is calculated using the following formula: In the formula: CESDR is the composite supply-demand ratio index. A positive CESDR indicates oversupply, a CESDR of 0 indicates supply-demand balance, and a negative CESDR indicates undersupply; S i D is a measure of the supply side of the i-th ecosystem service; i S is a measure of the demand side of the i-th ecosystem service; imax D represents the maximum measure of the supply side of the i-th ecosystem service; imax is the maximum measure of the demand side of the i-th ecosystem service; n is the total number of ecosystem services.

8. The method for identifying ecological protected areas based on the mismatch between supply and demand of ecosystem services according to claim 2, characterized in that, The driving factors include land use type, soil type, elevation, temperature, solar radiation, precipitation, potential evapotranspiration, vegetation index, vegetation available water content, GDP, and population density.

9. The method for identifying ecological protected areas based on ecosystem service supply-demand mismatch according to claim 2, characterized in that, The process also includes step S6: After identifying the ecological protection zone, if there is a conflict with the local economic development plan, the weights of various ecosystem services are adjusted, and then the weighted overlay tool of ArcGIS software is used to synthesize the data to generate a summary raster that reflects the spatial coordination / trade-off characteristics of ecosystem service supply and demand. Then, step S5 is repeated to complete the identification of the ecological protection zone.

Images