Identifying Critical Area of Ecosystem Service Supply and Demand at Different Scales Based on Spatial Heterogeneity Assessment and SOFM Neural Network

The purpose of this study was to establish a spatial multi-scale integrated assessment framework for critical areas of ecosystem service supply and demand, in order to provide theoretical support for regional ecological protection planning and re ﬁ ned management. Taking the typical hilly area of the upper reaches of the Yangtze River in China as an example, based on the assessment matrix of land use and ecosystem services, we used the method of spatial heterogeneity assessment and self-organizing feature mapping (SOFM) to explore the identi ﬁ cation and regionalization of critical areas of ecosystem services at regional and small scales. The results show that there was spatial heterogeneity and scale effects of ecosystem services under the two scales. The small-scale demand the importance of was follows: forest area in the upper part of the mountain > orchard and area in the middle and lower part of the mountain > valley farming area > ﬂ at town and farming area. The regional-scale supply – demand pattern was greatly affected by landscape structure, location conditions and social economy, the importance of services as follows: south > west > north > central. The SOFM network quantitatively identi ﬁ ed four types of ecological areas with clear dominant functions at regional and small scales. The balance between supply and demand in the studied ecosystem service areas was I < II < III < IV, in which IV was the critical area in terms of supply and I was the critical area in terms of demand. This assessment framework can improve the spatial accuracy and objectivity of the quanti ﬁ cation and mapping of ecosystem services, and provide new ideas for multi-scale identi ﬁ cation and expression of ecosystem services.


INTRODUCTION
Ecosystem services are products and services that human beings obtain directly or indirectly from ecosystem functions (Costanza et al., 1997;MA, 2005). The sustainable supply of ecosystem services forms the environment and the material basis for human survival and development, which has an extremely high or even incalculable value (Yan et al., 2017). However, continuous population growth, industrialization, and urbanization have placed considerable pressure on ecosystems, leading to the serious degradation of ecosystem service capacity. Over the past 50 years, provisioning services have increased globally, but most other types of services have shown a downward trend, and among these the declining rate of regulating services is extremely significant . On the other hand, the human socio-economic system has underestimated the value of ecosystem services for a long time, which has led to the excessive consumption of ecosystem services and serious imbalances and damage to ecosystems (Costanza et al., 1997;Fu and Zhang, 2014;Xie et al., 2015). The protection of natural ecosystems and enhancement of the quantity and quality of ecosystem services has become an urgent challenge globally (TEEB Foundations, 2010).
The relationship between supply and demand in relation to ecosystem services and their spatial measurement are the key to identifying regional eco-environmental problems and revealing their internal driving forces . The products and services supplied by ecosystems, together with human demand and consumption, constitute the dynamic process of ecosystem services, involving the natural ecosystem and the human socioeconomic system (Ma et al., 2017). The supply of ecosystem services is spatially constrained by biophysical and chemical processes, as well as the spatial distribution and structure of land use/cover (LUCC). The demand for ecosystem services is determined by the social and economic activities of individuals and groups, with regional differences and demand elasticity differences (Yan et al., 2017). In recent years, many scholars have conducted extensive research on the supply, demand, flow, and trade-offs of ecosystem services in different regions (Fisher et al., 2009;Burkhard et al., 2012Burkhard et al., , 2014Villamagna et al., 2013;Peng et al., 2017;Wang J. et al., 2019). The research objects were mainly concentrated on the supply and demand of a single ecosystem service type, with more research focusing on provisioning services, followed by regulation services and cultural services, and less attention to maintenance services (Wang Z. Z. et al., 2019). The research scale involves global, continental, national, regional and local scales, but the research has been mainly based on a single spatial scale (Ma et al., 2017), and there has been a lack of multi-scale research. Research cases are mainly concentrated in Europe, the United States and Asia (Wei et al., 2017). Among them, provisioning services analysis focuses on water and food Yang et al., 2012;Boithias et al., 2014;Jie et al., 2015); Regulation and maintenance services analysis focuses on flood regulation, pollination, and climate regulation (Stürck et al., 2014;Zhao and Sander, 2015); Cultural services analysis focuses on recreation (Beichler, 2015;Peña et al., 2015;Tratalos et al., 2016). On the other hand, land use change caused by human activities is the main factor that causes spatial differences and imbalances between supply and demand in relation to ecosystem services (Li et al., 2011). Burkhard et al. (2009Burkhard et al. ( , 2012 proposed the relationship matrix assessment method from the perspective of the internal relationship between land use and ecosystem services, which achieved the rapid assessment of the supply and demand potential of ecosystem services. This method can be used for ecosystem assessments with high complexity and uncertainty, and has attracted widespread attention and application (Li et al., 2016;Guan et al., 2020). However, the relationship matrix method cannot sufficiently reflect the internal heterogeneity and marginal effect of ecosystem service supply and demand (Ma et al., 2017), especially at small and micro-scales. Therefore, it is necessary to improve the relationship matrix method from the perspective of spatial heterogeneity to enhance its prediction accuracy in terms of the supply and demand distribution of ecosystem services at a small scale.
The identification of critical areas is the core of ecological protection zoning and planning, and the source of much difficulty in the practical application of ecosystem services theory. Previous studies often clustered and partitioned areas based on administrative or watershed units, but it is difficult to ensure the consistency of physical geography and ecological functional characteristics within these units, which seriously limits the application of multi-scale research (Mao et al., 2019). In recent years, some scholars have utilized hotspot analysis to explore identification methods for important areas in terms of multi-scale ecosystem services (Cai et al., 2017;Hou et al., 2018;Wang J. et al., 2019). However, although hotspot analysis method is suitable for one-dimensional or low-dimensional data, multi-dimensional ecosystem service data must be reduced and integrated before hotspot analysis, which increases the uncertainty of the results to a certain extent. Self-organizing feature mapping (SOFM) can classify multidimensional data sets based on their individual characteristics. It has the advantages of maintaining the topological structure, self-organization, self-adaptability, and strong fault tolerance. SOFM can be used to achieve bottom-up natural zoning, which is an effective quantitative method to connect the characteristics of physical geographical units and natural zoning (Ma et al., 2013). Some scholars have applied SOFM method to ecological function zoning (Ma et al., 2013;Mao et al., 2019) and ecological protection and restoration divisions (Wu et al., 2020;Xie et al., 2020), which has proved that the method is feasible for ecosystem services. However, these studies have primarily focused on the supply of ecosystem services at a single scale, rather than from the perspective of supply and demand. These applications are limited to identifying the critical areas of ecological protection, and there are still defects in the identification and characterization of critical areas relating to the balance of supply and demand in ecosystem services.
At present, China's per capita ecological assets are lower than the world average (Xie et al., 2017). In order to solve the environmental crisis, the Chinese government is trying to build an ecological protection decision-making system based on the "Ecological Redline" to ensure the integrity of critical ecosystems and the effective supply of important ecosystem services. The question of how to accurately identify the critical areas of ecological protection has become the primary problem to be solved. Recently, various studies have focused on Frontiers in Environmental Science | www.frontiersin.org July 2021 | Volume 9 | Article 714874 agriculture (Elahi et al., 2018a;Elahi et al., 2019aElahi et al., 2021a;Elahi et al., 2021b) and livestock sectors (Elahi et al., 2017;Elahi et al., 2018b;Elahi et al., 2019b) for estimation of resources and institutional barriers in the land use. Moreover, some studies focused on environmental issues Zhong et al., 2020); however, identification of the critical area of ecosystem service supply and demand at different scales has not been determined yet. Therefore, the main aim of this study is to propose a spatial multi-scale integrated ecosystem service pattern assessment framework to improve the accuracy of the prediction and identification of critical areas of the supply and demand of ecosystem services, and applies the framework to the ecological environment management of typical hilly areas in the upper reaches of the Yangtze River.

Classification of Land Use and Ecosystem Services
China's current land use classification system is based on the economic and social attributes of land, and is mainly used for land use management. This classification system does not reflect the ecological function of land very well. Specifically, it does not include ecological land types, and thus cannot meet the needs of ES expression, analysis and mapping (Chen and Shi, 2005;Deng et al., 2009;Zhang et al., 2017). Therefore, based on the identification and evaluation of land production, living and ecological functions, we established a new threelevel land use classification system in a previous study, which was named the production-living-ecological land (PLEL) classification system (Liao et al., 2019). This system included three first-level categories. Production land (PL) is related to the industrial structure and refers to the main function of providing industrial products, agricultural products, and service products. Living land (LL) is related to carrying and maintaining human settlements and refers to the main functions of providing human living, consumption, leisure, and entertainment. Ecological land (EL) is related to the natural ecosystem and refers to the main function of providing ecological products and services, as well as playing an important role in regulating, maintaining, and ensuring regional ecological security. Subsequently, 21 land use types were classified, according to the land cover and ecosystem types ( Table 1).
According to the different ecological functions, ecosystem services were divided into four categories: provisioning services, regulating services, habitat services, and cultural services (MA (Millennium Ecosystem Assessment), 2005; TEEB Foundations, 2010). Provisioning services are the material or energy outputs from ecosystems, such as food and fiber. Regulating services are the services that ecosystems provide by acting as regulators, such as climate regulation. Habitat services are the services that ecosystems provide as habitats for migratory species and for gene-pool protection, such as the maintenance of genetic diversity. Cultural services are the nonmaterial benefits that people obtain from ecosystems through spiritual enrichment, cognitive development, recreation, and aesthetic experiences.

Assessment Matrix of Supply and Demand of Ecosystem Services
Based on the non-monetary assessment method proposed by Burkhard et al. (2009Burkhard et al. ( , 2012Burkhard et al. ( , 2014 and the independent scores of regional experts, we established assessment matrixes of ecosystem services' supply and demand potentialities for different land use types (Figure 1, 2). In the assessment matrixes, 21 land use types are placed on the y-axis, and 11 ecosystem services are placed on the x-axis. The score of the matrix has the following meanings: 0 no relevant potentiality of supply or demand; 1 very low relevant potentiality of supply or demand; 2 low relevant potentiality of supply or demand; 3 medium relevant potentiality of supply or demand; 4 high relevant potentiality of supply or demand; and 5 very high relevant potentiality of supply or demand.

Spatial Heterogeneity Correction Model of Ecosystem Services
Due to the spatial heterogeneity of the influencing factors, such as natural, social, economic, and population factors, etc., the same land use types or ecosystems can have different supply and demand potential in terms of ES in different regions. Therefore, this study proposed an assessment method for the spatial heterogeneity of ES, which improves on the assessment matrix method of Burkhard et al. (2012).
1) In this study, biomass, soil, meteorology, topography, and other indicators with strong correlations or influences on ecosystem service supply were selected to establish the spatial heterogeneity correction model for each grid (10 m × 10 m) in the study area. The formula is as follows: where C ij refers to the spatial heterogeneity correction coefficient of ES supply potential i in grid j, including: food (C 1 ), materials (C 2 ), water (C 3 ), air quality regulation (C 4 ), climate regulation (C 5 ), water purification (C 6 ), regulation of water flows (C 7 ), erosion prevention (C 8 ), maintenance of soil fertility (C 9 ), habitat services (C 10 ), and cultural and amenity services (C 11 ); S ij refers to the spatial heterogeneity index of ecosystem supply service i in grid j; S i mean refers to the mean spatial heterogeneity index of ecosystem supply service i in the study area. To avoid the zero value of C ij , we assign C ij as being less than 0.1 to 0.1 after the initial calculation. The calculation method of the spatial heterogeneity index of 11 ES supply categories are shown in Table 2.
2) Demand for ecosystem services refers to the amount of ES consumed (available) or desired by human society. The demand of ES is linked to economic level, population density, and social preference. With reference to the relevant research literature Wang J. et al., 2019), three socio-economic indicators of city distance, population density, and gross domestic production (GDP) were selected to express the spatial heterogeneity of ES demand. City distance refers to the buffer value of cities at the county level and above, reflecting the intensity of human consumption of ES. Population density can reflect the quantity of demand for ES. The greater the population density, the greater the total demand for services. The GDP can indirectly reflect human preferences for ES. The higher the regional economic level, the higher the expectation for ES. The formula is as follows: Frontiers in Environmental Science | www.frontiersin.org July 2021 | Volume 9 | Article 714874 CDIS j cdis j cdis mean (3) PD j pd j pd mean (4) GDP j gdp j gdp mean (5) where X j refers to the spatial heterogeneity correction coefficient of ES demand potential in grid j; CDIS j refers to the normalization of the distance from grid j to the nearest city; cdis j refers to the distance (m) from grid j to the nearest city; cdis j mean refers to the mean distance (m) to the nearest city in the study area; PD j refers to the normalization of population density in grid j; pd j refers to the population density (persons/km 2 ) in grid j; pd mean refers to the mean population density (persons/km 2 ) in the study area; GDP i refers to the normalization of GDP in grid j; gdp j refers to the GDP (RMB 1000) in grid j; gdp mean refers to the mean GDP (RMB 1000) in the study area.

Supply and Demand Index Model of Ecosystem Services at the Small Scale
The purpose of small-scale identification is to reveal the spatial distribution characteristics of ecosystem service supply and demand inside the sample towns. Therefore, based on the assessment matrix and the spatial heterogeneity correction model of ecosystem service, the supply and demand index models were established to identify critical areas at the small scale. Then, we converted the PLEL vector maps of sample towns to raster maps, and calculated the supply, demand, and balance index of each grid. The model can be expressed as: where ESPI ij refers to the supply index of ecosystem service i in grid j; ESDI ij refers to the demand index of ecosystem service i in grid j; MSP ikj refers to the assessment matrix score of land use type k, corresponding to the supply potential of ecosystem service i ( Figure 2) in grid j; MSD ikj refers to the assessment matrix score of land use type k, corresponding to the demand potential of ecosystem service i ( Figure 3) in grid j; C ij , X j , and A j are the same as above; ESBI ij refers to the supply and demand balance index of ecosystem service i in grid j, ESBI > 0 means that supply exceeds demand, ESBI 0 means that supply and demand are in balance, and ESBI < 0 means that supply is less than demand. Frontiers in Environmental Science | www.frontiersin.org July 2021 | Volume 9 | Article 714874 5 TABLE 2 | Spatial heterogeneity index model of the supply of 11 ecosystem services.

Ecosystem service types Model Description
Food (S 1j ) S 1j PCROP j There is a strong correlation between food supply capacity and land productivity potential (Wang et al., 2015). The index of food supply capacity is estimated by using the spatial distribution data of farmland productivity potential in 2015. PCROP j refers To value of farmland productivity potential in grid j (kg/ha) Materials (S 2j ) S 2j NDVI j Based on the correlation between the vegetation index and timber volume (Wang et al., 2007), the normalized difference vegetation index (NDVI) is used to represent the spatial heterogeneity of raw materials. NDVI j refers To value of annual NDVI in grid j Water (S 3j ) S 3j (1 − AET j /P j ) × P j The spatial heterogeneity of water supply is expressed by the water yield, which is assessed based on (Budyko, 1974) water-heat coupling equilibrium hypothesis and annual average precipitation data. AET j refers To the value of annual actual evapotranspiration in grid j (mm); P j refers to value of average annual precipitation in grid j (mm) Air quality regulation (S 4j ) S 4j LAI j Vegetation can significantly reduce atmospheric gaseous and solid pollutants, and its ability to absorb pollutants is related to leaf area and structure (de Groot et al. 2010;Yi et al., 2017;Xu et al., 2018). Therefore, we use the leaf area index (LAI) to represent the spatial heterogeneity of air quality regulations. LAI j refers To the value of annual average LAI in grid j Climate regulation (S 5j ) S 5j G veg j + G soil j Carbon pools in terrestrial ecosystems play an important role in controlling global climate change (Yu et al., 2010). According to the methods of Specifications for the Assessment of Forest Ecosystem Services (Chinese standard LY/T 1721-2008), the spatial heterogeneity of climate regulation is expressed by carbon fixation. G veg refers To carbon fixation by vegetation in grid j (t/year); R is the carbon content in CO2, which is 27.27%; A j refers to the area of grid j (ha); B j refers to the value of annual net primary productivity (NPP) in grid j (t/ha/year); G soil refers to carbon fixation by soil in grid j (t/year); H refers to the thickness of surface soil (20 cm is used in this study) (cm); OM j refers to the content of soil organic matter in grid j (%); W j refers to the soil bulk density in grid j, which is 1.3 according to the average value in the study region (g/cm3) Relevant studies have shown that water network density has a significant positive impact on water purification services . Therefore, we use water network density to represent the spatial heterogeneity of water purification services. WND j refers To the value of water network density in grid j (ha/km2) Regulation of water flows (S 7j ) According to the NPP quantitative index method of the Technical Guidelines for the Delineation of Red Lines for Ecological Protection (2015) issued by the Ministry of Environmental Protection in China, we use the water conservation assessment method to evaluate the spatial heterogeneity of the regulation of water flows. B j and P j is the same as above; FSIC j refers to value of the soil penetrating capacity factor in grid j; FSIC max and FSIC min refers to the maximum and minimum values of the soil penetrating capacity factor in the hilly area of the study region; P max and P min refers to the maximum and minimum values of average annual precipitation in the hilly area of the study region; Eq. 12 is the empirical formula of the soil penetrating capacity factor (Che 1995); SAND j refers to the soil sand content in grid j (%); SLOP j refers to the slope in grid j (degree); SLOP max and SLOP min refer to the maximum and minimum values of the slope in the study region According to the NPP quantitative index method of the Technical Guidelines for the Delineation of Red Lines for Ecological Protection (2015) issued by the Ministry of Environmental Protection in China, we use the erosion prevention assessment method to evaluate the spatial heterogeneity of erosion prevention. B j , SLOP j , OM j and SAND j are the same as above; K j refers to value of the soil erosion factor in grid j; SILT j refers to the soil silt content in grid j (%); CLAY j refers to the soil clay content in grid j (%) The spatial heterogeneity of the maintenance of soil fertility is expressed by the comprehensive index of soil nitrogen, phosphorus, and potassium contents (Liu et al., 2018). SN j , SP j , and SK j refer to the soil's available nitrogen, phosphorus. And potassium content in grid j (mg/kg) Habitat services (S 10j ) According to the NPP quantitative index method of the Technical Guidelines for the Delineation of Red Lines for Ecological Protection (2015) issued by the Ministry of Environmental Protection in China, we use the biodiversity assessment method to evaluate the spatial heterogeneity of habitat services.
(Continued on following page) Frontiers in Environmental Science | www.frontiersin.org July 2021 | Volume 9 | Article 714874 6 Supply and Demand Index Model of Ecosystem Services at the Regional Scale The purpose of regional-scale identification is to reveal the critical areas of the supply and demand of ecosystem services in the study area, with sample towns as the basic unit. Therefore, based on the supply and demand index models at the small scale, the total supply and demand index models were established to identify critical areas at the regional scale. Due to the different areas of the sample towns, we used the unit area index to revise the supply, demand, and balance index of each town. The model can be expressed as: where ESBT in ′ refers to the total supply and demand balance index of ecosystem service i of sample town n, m is the number of grids in sample town n, and A n refers to the area (ha) of sample town n. ESBT n refers to the total supply and demand balance index of sample town n, ESBT > 0 means that supply exceeds demand, ESBT 0 means that supply and demand are in balance, and ESBT < 0 means that supply is less than demand.

Identifying Critical Areas of Ecosystem Service Supply and Demand
Self-organizing feature mapping (SOFM) is an unsupervised artificial neural network, which can map a high-dimensional data set to a low-dimensional space and obtain the similarity relationships between data. The SOFM network consists of two parts: the input layer and the output (competition) layer ( Figure 3). The dimension of the input vector is the same as that of each node in the competition layer (Kohonen, 1997). Therefore, data clustering results can be obtained on the premise of keeping the topological structure of data set unchanged, and this is an effective method to identify and regionalize ecosystem services (Ma et al., 2013). The learning algorithm process of the SOFM network is as follows: 1) Initialization, initialize the connection weights from N input neurons to output neurons, assign a small random initial value w ij , and set an initial field for any output node j. 2) Provide a new input vector X. 3) Calculate the Euclidean distance d j between the input sample and each output neuron j. Then, select the neuron k with the smallest distance as the winning output unit. The calculation formula is: 4) Given a surrounding domain S k (t) to update the weight of the winning node and the nodes in the surrounding domain, the amount of weight change is: where, η refers to positive learning rate, decreasing with time.

Ecosystem service types
Model Description B j , P j , P max , and P min are the same as above; TEM j refers to value of ghd average annual temperature in grid j; TEM max and TEM min refer to the maximum and minimum values of average annual temperature in the hilly area of the study region; ALT j refers to value of altitude in grid j; ALT max and ALT min refers to the maximum and minimum values of the altitude in the study region Cultural and amenity services (S 11j ) The spatial heterogeneity of cultural and amenity services is expressed by the kernel density value of tourist attractions (Ouyang et al., 2013), which are generated by means of kernel density analysis using ArcGIS 10.6. r refers to the default search radius of the grid j (m); n refers to the number of tourist attractions within the search radius distance of the grid j; d j is the distance between grid j and a certain tourist attraction (m) Frontiers in Environmental Science | www.frontiersin.org July 2021 | Volume 9 | Article 714874 5) Provide a new input vector and repeat the above learning process.

Study Area
The hilly area of the Sichuan Basin is located in the upper reaches of the Yangtze River in Southwest China (28°10′-31°53′ N, 102°47′-107°28′ E) with an area of about 84,000 km 2 , including 70 counties (Figure 4), which is the most typical hilly area of the Mesozoic continental red strata in China. This region is an important functional region of urban development, food production and ecological protection in Western China and the upper reaches of Yangtze River, which plays an important role in the overall situation of regional economic and social development. However, due to the fragmentation of the natural landscape, low forest coverage, and low efficiency of agricultural production and land use, serious problems such as soil degradation, environmental pollution, resource depletion, and a decline in biodiversity have arisen, with this area becoming the most serious area of soil erosion in the upper reaches of the Yangtze River (Qiao and Yan, 2008;Wei et al., 2009). Because it is difficult to collect the land use/cover and socio-economic data for all towns of the study area, we therefore selected 215 typical hilly towns in the north, middle, west, and south of the Sichuan Basin as samples for this case analysis.

Data Sources and Processing
In this study, ERDAS IMAGINE software was utilized to enhance and correct the Google Earth satellite images. Referring to the production-living-ecological land (PLEL) classification system, the PLEL vector database (1:10,000 scale) of sample towns in the study area was established through field investigations and the visual interpretation of the remote sensing images using ArcGIS 10.6. Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTERGDEM) data were gathered from NASA. The data set of farmland productivity potential, the normalized difference vegetation index (NDVI), net primary productivity (NPP), soil texture, annual average precipitation, annual average temperature, distribution of tourist attractions, and GDP in 2015 were provided by the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC) (http://www.resdc.cn).
Frontiers in Environmental Science | www.frontiersin.org July 2021 | Volume 9 | Article 714874 9 organic matter, alkali-hydro nitrogen, rapidly-available phosphorus, and rapidly-available potassium was calculated through Kriging interpolation of 38,963 soil sampling points in the soil testing and formula fertilization database of Sichuan Province. Socioeconomic data were collected from the Statistical Yearbook of Sichuan Province and the Statistical Yearbook of various counties (2015).

Characteristics of Ecosystem Service Supply and Demand
Based on the spatial heterogeneity assessment of ecosystem services, 11 supply and demand balance indices of ecosystem services (ESBI values) were calculated. Figure 5 shows that food (E1) supply was generally higher than demand in the study area, which was the main agricultural production area in Sichuan Basin, with more than 65% cultivated land. Spatially, the high value of food services decreased gradually from the low hilly area in the middle to the surrounding high hilly area. Correspondingly, high population density and low vegetation coverage led to a significant reduction in the supply capacity of ecosystem services in the central and northern parts of the study area. Therefore, the ESBI value of E2-E10 ecosystem services from high to low was as follows: south > west > north > middle. In addition, cultural and amenity services (E11) were greatly affected by socio-economic and location conditions, so the supply capacity of ecological areas closer to cities or settlements was generally higher. Figure 6 shows that the supply and demand of ecosystem services showed obvious spatial heterogeneity and regularity at different scales in typical hilly areas of SW China. At the regional scale ( Figure 6, right), the total supply potentialities of ecosystem services increased gradually from the northeast to the southwest, which was consistent with the distribution trend of natural environmental factors, such as better vegetation coverage and hydrothermal conditions. On the contrary, the total demand potentialities of ecosystem services decreased from the northeast to the southwest. Especially in the middle and northeast of the study area, due to the high population density and the agriculture-dominated rural areas, the demand for ecosystem services was much higher than the supply. At the small scale, we take Qiuchang town, in the south of the study area, as an example for comparative analysis (Figure 6, left). The ESBI map, based on the spatial heterogeneity assessment, clearly expressed the spatial distribution of the critical area of ecosystem service supply and demand at the small scale. The ESBI of Qiuchang town was mainly affected by topography and land use. The supply area of ecosystem services is mainly distributed in the central and southeast mountainous areas, whereas the demand area of ecosystem services is mainly concentrated in the agricultural and urban areas of the valley. Figure 7 shows that the supply and demand potential of ecosystem services at different scales was unbalanced, and the structure was similar. The supply of food, materials, air quality regulation, climate regulation, habitat and cultural and amenity services was greater than the demand, whereas the demand for water, water purification, erosion prevention and maintenance of soil fertility was significantly higher than the supply. At the regional scale ( Figure 7A), the total supply met nearly 91% of the total demand. The supply potentialities of ES were mainly composed of regulating and provisioning services, at 42.66 and 36.13% of the total supply, respectively, which was mainly provided by arbor forest land, riverine wetlands, paddy fields, and orchards. Meanwhile, the demand potentialities of ES were Frontiers in Environmental Science | www.frontiersin.org July 2021 | Volume 9 | Article 714874 also composed of regulating and provisioning services, at 77.57 and 19.41% of the total demand, respectively. At the small scale ( Figure 7B), the total supply of Qiuchang town met 175% of the total demand. The supply potentialities of ES were also mainly composed of regulating and provisioning services, which accounted for 52.94 and 27.39%, respectively, and the demand potentialities of regulating and provisioning services accounted for 77.96 and 19.12% respectively.

Identification and Regionalization of Ecosystem Service Supply and Demand
In this paper, the SPSS modeler was used to cluster the SOFM network to identify the critical areas of the supply and demand of ecosystem services. At the small scale, the ESBI values of 11 ecosystem services of each grid were used as the input parameters, the training samples were 893,887 grid units of Qiuchang Town; at the regional scale, the ESBT' values of 11 ecosystem services of sample towns were used as the input parameters, and the training samples were 215 typical towns in hilly areas; the input parameters were trained for 1,000 times to build a stable SOFM model; the other parameters of the model used the software default values, and the output nodes were 12. In order to maintain the differences between different types, combined with the characteristics of regional ecosystem services, the SOFM clustering results are combined into four types. The calculation results of SPSS modeler software showed that the condensed and separated contour values at the two scales were both greater than 0.5, and the classification results were reliable. In addition, the spatial weight order of the small scale was: E9 > E8> E7 > E5 > E4 > E10 > E3 > E2 > E1 > E6 > E11; the spatial weight order of the regional scale was: In general, regulation and provision of services had a greater impact on ecosystem services in the study area. According to the clustering results of the SOFM network, the critical areas of ecosystem service supply and demand were identified through four categories at two scales. Figure 8 shows the spatial distribution of supply and demand types of ecosystem services at the regional scale (right) and the small scale (left). Meanwhile, we analyzed the characteristics of ecosystem services in different regions, and Figure 9 shows the average value of the supply and demand balance indices of 11 ecosystem services.

1) Regional Scale
Type I: This type accounted for 32.09% of the sample towns, mainly distributed in the middle and northeast of the study area. This type of area was relatively flat land, close to the city, with high population density, a high proportion of agricultural and residential land, but the area of forest and wetlands was small, and the landscape was fragmented. Therefore, the imbalance between supply and demand was the most serious in this region, and it showed a strong demand for most ecosystem services, except for food and cultural and amenity services. Due to the high intensity of human disturbance in the region, it is necessary to strictly control the utilization of natural resources, implement soil and water conservation projects, and strengthen agricultural ecological protection and vegetation restoration in these areas.
Type II: This type accounted for 30.70% of the sample towns, mainly distributed in the transition zone around type I. The ecosystem service structure of this region was similar to that of type I, except that the supply and demand of habitat and materials services were balanced. The natural erosion of the surface was also serious, so it is necessary to prevent unreasonable land use from aggravating soil erosion, such as through steep slope reclamation.
Type III: This type accounted for 23.72% of the sample towns, mainly distributed in the west and south of the study area. The supply capacity of ecosystem services in this region was higher than the demand, and only water, water purification, and the maintenance of soil fertility were unable to meet the demand. It is necessary to strengthen the protection of natural vegetation with the function of water conservation and expand the conversion of farmland to forest and grasslands in these areas.
Type IV: This type accounted for 13.49% of the sample towns, mainly distributed in the south and southwest of the study area, with good natural environment quality and strong ecosystem service supply capacity. The area of forest and wetlands was large, and the landscape was concentrated and continuous. It is the critical area of ecosystem service supply at the regional scale, and Frontiers in Environmental Science | www.frontiersin.org plays an important role in the optimal management of ecosystems. In the future, the region should strengthen the protection of natural ecosystems and biodiversity, and should be included in the regional ecological red line for key management.

2) Small Scale
Type I: Areas of this type accounted for 23.72% of Qiuchang Town, mainly distributed in the urban and rural residential areas, with strong demand for ecosystem services. The area mainly provided industrial, residential, cultural, and entertainment services. The supply of ecosystem services cannot meet the needs of the people, and it was the sink of the ecosystem service demand in Qiuchang. Therefore, it is necessary to strengthen ecological restoration and reconstruction, optimize land use structure, and construct ecological network systems such as corridor and buffer zones, to ensure the effective flow and supply of ecosystem services.
Type II: Areas of this type accounted for 39.39% of Qiuchang Town, mainly distributed in gentle valleys and lower hills. The land use was mainly cultivated land, and the surface cover was crops, which was the critical area of food supply in Qiuchang. Single and unbalanced ecosystem services were the obvious FIGURE 8 | Regionalization of ecosystem service supply and demand at different scales in typical hilly areas, SW China. Left: small scale; right: regional scale.
Frontiers in Environmental Science | www.frontiersin.org July 2021 | Volume 9 | Article 714874 characteristics of this region. The geographical conditions and hydrothermal resources of this region were good, and the potential for land development and utilization was great. However, it is still necessary to reasonably organize agricultural production and rural economic activities, strengthen the management and protection of agricultural ecosystems, prevent agricultural non-point source pollution, and enhance the area's ability to resist natural disasters.
Type III: Areas of this type accounted for 12.00% of Qiuchang Town, mainly distributed in the transition zone between mountains and valleys. The land use was mainly orchard and arid land in the middle of hills. The supply and demand of this region were generally balanced, and the supply of 10 ecosystem services (except water) exceeded demand. This region was mainly a semi-natural ecosystem with a high altitude and large slope, and showed great potential to improve ecosystem services. Through clean production, slope farmland management, and conversion of farmland to forest, improving the ecological functions of soil and water conservation will be the key issues of ecosystem management.
Type IV: Areas of this type accounted for 42.42%, mainly distributed in mountainous areas with high forest coverage, which was the critical area of ecosystem services supply in Qiuchang town. The good natural environment and biodiversity of this region provide an important basis for ecological protection. It is necessary to strictly control the development and construction within and around this region, to reduce the human disturbance of the ecosystem, and finally to establish this area as a key zone of ecosystem conservation.

DISCUSSION
Identifying critical areas of ecosystem service supply and demand is very important for ecological protection planning, resource allocation, and management decision-making. The primary problem in critical area identification is the quantitative assessment of ecosystem service supply and demand. Due to the limitation of small-scale applications of Burkhard's method, this paper constructed a spatial heterogeneity assessment model, and attempted to accurately simulate the supply and demand of ecosystem services based on grid units. The case study results show that the grid-based spatial heterogeneity correction method can reduce the uncertainty of the assessment matrix method to a certain extent (Ou et al., 2018), and provide a reference for the accurate identification and expression of ecosystem services. On the other hand, due to the subjectivity of regions and experts, the supply and demand scores of different regions had certain differences. However, from the construction of the evaluation matrix, the scoring trends of this paper and related research were generally consistent (Cai et al., 2017;Ou et al., 2018), and the research results were still comparable. In addition, the study on the scale response characteristics of ecosystem services (Huang et al., 2019) showed that the township scale was the characteristic variation scale of ecosystem services, which can be used as an ideal scale level for ecosystem service evaluation and optimization. This result was consistent with our conclusion, indicating the rationality of choosing the township scale as a small-scale area study. However, the grid sizes currently selected for the study of ecosystem services at the township scale were quite different, so the scale issue needs to be further explored.
In this paper, a supply and demand index model of ecosystem services was constructed, and the unified quantification and mapping of ecosystem services based on multi-scale fusion were explored. In the case study, the small-scale supply-demand pattern was found to be affected by microtopography and land use patterns, whereas the regional scale supply-demand pattern was found to be affected by landscape structure, location conditions, and social economy. These results were consistent with the related research results on the regional scale (Yang et al., 2018;Xie et al., 2020) and the township scale (Jiang et al., 2020), indicating that the research method in this paper can objectively express and reflect the distribution characteristics of the supply and demand of ecosystem services at different scales. However, production and benefit areas of ecosystem services are not completely coincident in space (Fisher et al., 2009); our method cannot effectively express the flow characteristics of ecosystem services. Accurate simulation of the spatial impact range and transmission path of various ecosystem services will be the focus of future research. The SOFM network model was used to identify and extract the critical areas of ecosystem service supply and demand, while avoiding the subjectivity of traditional methods. The case study shows that the identification results of critical areas at different scales showed clear characteristics of dominant ecological function and the supply-demand relationship. Meanwhile, the differences in the functional characteristics of each type of region at the small scale was more significant than that at the regional scale. This may be due to the fact that there were more training samples at the small scale, so the constructed SOFM model had higher accuracy and better recognition ability at this scale. However, SOFM neural network clustering with multidimensional raster data often has some problems, such as high fragmentation of the edge region and unclear boundaries (Ma et al., 2013;Mao et al., 2019), which limits the practicability of the research results. Effectively identifying and merging the clustered, broken, and mixed blocks into complete blocks remains an urgent problem.

CONCLUSION
This paper presents a new framework for the identification and expression of the supply and demand patterns of ecosystem services. Through this framework, we identified the critical areas of ecosystem service supply and demand in the typical hilly area of the upper reaches of the Yangtze River, and provide a reference for the delimitation of ecological protection red lines and ecosystem management decisions at different scales.
First, a quantitative method of assessing spatial heterogeneity was applied to reveal the spatial variation characteristics of the supply and demand of ecosystem services, which extended the theoretical system and application scope of Burkhard's method, and improved the accuracy of the identification and expression of Frontiers in Environmental Science | www.frontiersin.org July 2021 | Volume 9 | Article 714874 ecosystem services. Case studies show that there was spatial heterogeneity and scale effects of ecosystem services under the two scales. The small-scale supply-demand pattern was greatly affected by microtopography and land use patterns, and the importance of ecosystem services was as follows: forest area in the upper part of the mountain > orchard and dry area in the middle and lower part of the mountain > valley farming area > flat town and farming area. The regional-scale supply-demand pattern was greatly affected by landscape structure, location conditions and social economy, and the importance of ecosystem services was as follows: south > west > north > central. Second, based on the SOFM network model, the identification and extraction of critical areas of the supply and demand of ecosystem services can avoid the subjectivity of traditional methods. The balance between supply and demand in the studied ecosystem service areas was I < II < III < IV, in which IV was the critical area in terms of supply and I was the critical area in terms of demand. The ecological function characteristics of each type of zoning were obvious, the spatial distribution characteristics were consistent with the reality, and the results were reliable.
Third, eco-environmental management policy recommendations: Type I areas need to implement a green, low-carbon economy regional development strategy, strengthen ecological restoration and reconstruction, and control the excessive consumption of natural resources. Type II areas need to implement a regional development strategy for cleaner production, optimize agricultural production and rural economic activities, and improve the ability to withstand natural disasters. Type III areas need to implement a regional ecological protection strategy of returning farmland to forests and grasslands to further improve the quantity and quality of natural ecosystems. Type IV areas need to implement an ecological red line management strategy, protect the integrity of the natural ecosystem and biodiversity, strictly control the development and construction of the interior and surrounding areas, and reduce human interference.
Due to the fragmentation and blurring of the boundaries of some critical areas of ecosystem services at small scales, effectively identifying and merging broken and mixed blocks into complete and clear partition boundaries will be the focus of future research. Although this research framework still has some limitations and uncertainties, it was able to realize the unified quantification, identification, and mapping of critical areas of ecosystem services at multiple scales and provide valuable information for the precise management of the ecological environment.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
Guitang Liao, Xuesong Gao and Liangji Deng conceived and designed the experiment. Zhengyu Lin and Peng He collected the data and analyzed the data. Wei Zhou and Chenghua Xu helped the language correction; Guitang Liao wrote the paper.

FUNDING
This work is financially supported by the Chengdu University of Information Technology (KYTZ202116).