Measuring the Relationship between Physical Geographic Features and the Constraints on Ecosystem Services from Urbanization Development

: Exploring the constraint relationship between physical geographic features and urbanization on ecosystem services is important for managing and optimizing regional ecosystem services. Taking Anhui Province as an example, we assessed the spatial and temporal evolution characteristics of ﬁve types of ecosystem services (habitat support, water production services, soil conservation, NPP, and carbon ﬁxation) and ﬁve types of urbanization levels (population, economic, social, ecological, and spatial) in 2000, 2010, and 2020, and integrated the constraint line method, bivariate spatial autocorrelation model, and spatial regression model to measure the relationship between ecosystem services. The spatial constraints between ecosystem services and urbanization level and natural topography in Anhui Province were measured using the constraint line method, bivariate spatial autocorrelation model and the spatial regression model. The results show that: (1) the spatial distribution of the ﬁve types of ecosystem services in Anhui Province is characterized as “low in the north and high in the south”. At the provincial level, the ﬁve ecosystem services in southern and central Anhui Province are synergistic, while the ﬁve ecosystem services in northern Anhui Province show a trade-off; (2) topography has different effects on the ﬁve ecosystem services with “exponential” effects on water production services and NPP, “positive convex” effects on habitat support, and “positive convex” effects on habitat support”; (3) the bivariate global autocorrelation Moran’s I index between ecosystem services and urbanization level in Anhui Province is signiﬁcant, conﬁrming that ecosystem services and urbanization are spatially related, where the development of population urbanization, spatial urbanization, economic urbanization, and social urbanization leads to the decrease in ecosystem services, and ecological urbanization promotes the increase in ecosystem services. In the spatial regression model, the Spatial Lag Model passed the signiﬁcance test, indicating that there is a spatial spillover effect between ecosystem services and urbanization. That is, changes in ecosystem services are inﬂuenced not only by their own urbanization elements, but also by urbanization elements in neighboring units or more distant units. Exploring the constraints of ecosystem services and identifying their interaction with urbanization can provide a scientiﬁc basis for land-use optimization, adjusting management measures and achieving regional sustainability.


Introduction
Ecosystem services contain numerous benefits, such as natural environmental conditions and utilities, that humans derive, directly or indirectly, from ecosystems and their ecological processes, and provide the resources and environmental basis for human survival [1,2]. However, in the context of rapid urbanization, uncontrolled human claims on natural resources and unscientific resource use have degraded nearly 60% of global terrestrial ecosystem services [3]. In 1987, the United Nations released the study "Our Common Future", which introduced the concept of sustainable development [4]. Economic growth, social progress and environmental protection were clearly defined as the three pillars of sustainable development. By the end of 2020, 56.1% of the world's population lived in cities. Cities occupy 3% of the world's land, consume 75% of the world's natural resources, emit 50% of the world's waste and 70% of the greenhouse gases, and also generate more than 80% of the GDP. Cities have become a key factor in advancing sustainable human development. Therefore, a scientific understanding of the interaction between ecosystem services and urban development and an in-depth investigation of its "mechanism-process" mechanism of influence have become an important means of theoretical support for promoting sustainable social development.
In recent years, scholars at home and abroad have discussed the spatial and temporal evolution of ecosystem services [5][6][7], influencing factors [8,9] and trade-offs and synergistic relationship measurements [10][11][12]. A large number of studies have been conducted in recent years. For example, Geng Tengwei et al. used geographic probes to find that urbanization is a key factor affecting the value of ecosystem services, and used geographically weighted regression models (GWR) to explore the spatially divergent characteristics of the intensity of urbanization effects [13], and Zhang Yushuo et al. [14] used multiple regression models to quantitatively reveal the influence of three dimensions, land use, social and economic dimensions, on the trade-offs and synergies of ecosystem services. Li Cheng et al. [15] presented a study that used the root-mean-square error index and redundancy analysis to investigate the factors influencing ecosystem services in the Yangtze River Delta urban agglomeration, and found that the increase in construction land was an important factor affecting the decline in ecosystem services; Howe et al. [16], based on existing studies on the trade-offs and synergistic relationships of ecosystem services, found that trade-offs in ecosystem services tend to be static, while natural environmental changes, ecosystem feedback and food webs may lead to lags in shifting trade-offs in ecosystem services. Alison et al. [17] assessed the complex relationship between ecosystem services and agricultural production and found that agriculture can be a source of degradation for many ecosystem services, such as loss of wildlife habitat, nutrient runoff, watercourse deposition, greenhouse gas emissions, and pesticide poisoning by humans and non-target species. Natural geographic features have been used by researchers to explore the spatial variations in ecosystem services as important influencing factors, e.g., Yang Sun et al. [18] assessed the value of ecosystem services in the middle reaches of the Yangtze River as a logarithmic function of topographic gradient, and the value of ecosystem services was greater in mountainous areas than in plains. The authors of [19] took the western plains and eastern mountains of the South Four Lakes Basin in eastern China as an example, and proved that topography is an important factor affecting the spatial distribution of multiple hydrological ecosystem services through correlation analysis and redundancy analysis. Researchers have studied the impact factor measures of ecosystem services such as GWR and geographic probes, and used a redundancy analysis to directly express the causality and intensity of the impact, i.e., determine the heterogeneity of positive and negative impact relationships or the spatial impact magnitude of a factor on ecosystem services. Numerous studies have provided empirical evidence, showing [20][21][22][23] that the relationship between ecosystem services and influencing factors changes with spatial scale and temporal scale, and shows nonlinear influence characteristics. The degree of influence that these factors have on ecosystem services in nonlinear processes (types of constraint curves, thresholds and eigenvalues) plays an important role in managing and balancing the negative impacts of ecosystem services and human activities [24][25][26]. Among the many factors affecting changes in ecosystem services, physical geographic features and human construction activities are considered to be the key driving forces that change and influence the spatial pattern of ecosystem services and trade-off relationships, meaning that physical geographic features and urbanization development are a long-standing focus of research on ecosystem service changes [27,28]. Theoretically, measuring the constraints on ecosystem services regarding physical and geographic features and urbanization development is important for enhancing the understanding of the drivers and mechanisms of ecosystem services [29,30]. To sum up, scholars pay more attention to the traditional linear regression and correlation analysis methods to measure the relationship between variables, but ignore the limiting effect of limiting variables on response variables in complex ecosystems affected by multiple factors. Scientific measurement of the limitations of natural geographical characteristics on ecosystem services is conducive to accurately understanding the marginal characteristics of the impact of urbanization on ecosystem services.
Anhui Province is located in the Yangtze River Delta city cluster, one of the six largest city clusters in the world, which is an important intersection between the Belt and Road and the Yangtze River Economic Belt, and an important platform for China's participation in international competition. The urbanization rate in Anhui Province has increased from 18.55% in 2000 to 58.33% in 2020, and the rapid urbanization activities have caused serious habitat fragmentation within the region, with an increase of 1.65% in construction land, a decrease of 0.5% in forest and grass land, and a decrease of 2.5% in arable land within the last 20 years. The spatial pattern of land cover in Anhui Province significantly varies with the changes in topography, and the altitude of the whole province gradually decreases from south to north. The southern part is mostly forested and grassland, while the northern part is mostly arable land and construction land, so it shares commonality and characteristics with cities in the lower Yangtze River basin. Against the background of the research presented above, this paper plans to conduct research focusing on three aspects: (1) At the research scale, this paper selects county administrative units (excluding municipal districts) to obtain measurements and statistics on ecosystem services and urbanization, mainly because the long-standing spatial raster data on population and economy are limited to the research scale, with limited spatial and temporal accuracy, an inconsistent data structure, and a mismatch of spatial units. (2) This paper aims to conduct a comprehensive evaluation of ecosystem services in Anhui Province during 2000, 2010 and 2020, based on the InVEST model and CASA model, and to evaluate the development level of urbanization in counties in Anhui Province during 2000, 2010 and 2020 in five dimensions: population, economy, society, space and ecology. (3) This paper aims to explore the constraints that natural topography places on ecosystems and the constraints threshold with the help of the constraint line method, based on the bivariate spatial autocorrelation model. The spatial interactions and spatial dependence between urbanization and ecosystem services in the study area were investigated based on Bivariate Moran's I, Ordinary Least Square (OLS), Spatial Lag Model (SLM), and Spatial Error Model (SEM), and the spatial interaction between urbanization and ecosystem services in the study area.

Overview of the Study Area
Anhui Province, located in East China, is a sub-center city of the Yangtze River Delta (Figure 1), between 114 • 54 -119 • 37 east longitude and 29 • 41 -34 • 38 north latitude, with a total area of 140,100,000 km 2 . There are 105 counties and districts, including 60 counties and 55 municipal districts. The topography of Anhui Province is diverse, with the Yangtze River and the Huai River crossing the whole territory from west to east, dividing the province into three natural regions: HuaiBei, Jianghuai and Jiangnan. The traditional zoning concept is divided into northern Anhui, central Anhui and southern Anhui. North Anhui refers to the area north of the Huai River, including Suizhou, Huabei, Bengbu, Fuyang, Huainan, Bozhou, with six cities; Central Anhui refers to the area north of the Yangtze River and south of the Huai River, including Hefei, Liuan, Chuzhou, Anqing, with four cities. South Anhui refers to the area south of the Yangtze River, including Huangshan, Wuhu, Maanshan, Tongling, Xuancheng, Chizhou, with six cities. To the north of the Huai River, the terrain is vast and open, and is part of the Great Plains of North China. The central part lies between Jiang and Huai, with mountainous hills. The banks of Yangtze River and around Chaohu Lake are low and flat, belonging to the famous middle and lower reaches of the Yangtze River Plain; the south is dominated by mountains and hills. It is located in the transition area between the warm temperate zone and subtropical zone. North of the Huai River is a warm, temperate, semi-humid monsoon climate, and south of the Huai River is a subtropical, humid monsoon climate, compatible with the north and south.

Data Sources and Processing
A rainfall dataset and evapotranspiration dataset from the National Earth System Science Data Center (http://www.geodata.cn/, accessed on 10 October 2021) were used to input rainfall and evapotranspiration data into the water production module of the InVEST model. The Digital Elevation Model (DEM) was obtained from the Geospatial Data Cloud website (http://www.gscloud.cn/, accessed on 20 July 2021) with a raster size of 90 m. Soil data were obtained from the Cold Region Dry Zone Science Data Center (http://data.casnw.net/portal/, accessed on 5 May 2021). Soil data were obtained from the Chinese Soil Dataset (v1.1), which was used to calculate the soil erodibility factor (K) in the soil conservation services, as well as the plant-available water content (PAWC) and depth-to-root restrictions in the water production module of the InVEST model [31]. The data of urbanization indicators were obtained from the Statistical Yearbook of Anhui Province, county statistical yearbooks and bulletins. The entropy value method was used to calculate the weights of each indicator, and the entropy value method was calculated by SPSSAU (https://www.spssau.com/, accessed on 5 May 2021). This was interpolated to 100 m and the coordinates were unified as WGS_1984.

Research Ideas and Methods
The InVEST model [32] was used to calculate habitat support services (HQ), soil conservation services (SDR), water production services (WY), and carbon fixation (CF) in the study area in 2000, 2010, and 2020. The CASA model was also used [33]. The NPP of the study area in 2000, 2010, and 2020 was calculated to replace the organic matter production function, and the five types of ecosystem services were weighted and calculated until the integrated ecosystem services. (2) Based on the statistical yearbooks and statistical bulletins of the study area in 2000, 2010 and 2020, we obtained the relevant evaluation indicators, and used the entropy value method to calculate the indicators at all levels to obtain the spatial and temporal evolution patterns of the five types of urbanization levels and the comprehensive urbanization level in the study area. (3) Based on the constraint line method [34,35], the spatial regression model was used to measure the spatial constraint relationship between integrated ecosystem services and five urbanization spatial subtypes( Figure 2).

1.
Habitat Quality (Habitant Quality) The InVEST model is used for ecosystem service function assessment, in which the habitat quality evaluation module is based on the link between land cover and habitat threat sources. It calculates the intensity of threat sources by considering the radius of stress, spatial weights, spatial attenuation types, etc., and combining the habitat adaptation of other land types and sensitivity to threat sources to obtain the habitat quality of the area using the following equation.
where Q xj represents the type of landscape in the study area, while j in x represents the habitat quality index of the raster H j . The value range of [0, 1] represents the habitat suitability score of the landscape type j in the habitat suitability scores; k is the halfsaturation constant, which is set according to the accuracy of data in the study area and is 50.5 in this paper; z. In this study, the Habitat Quality module parameter table was set based on the InVEST model manual and related studies.

Water Yield (WY)
In this study, the water production service refers to the surface water production in a certain area. The calculation of surface water production is based on a simplified hydrological cycle model that ignores the influence of groundwater and is determined by a combination of numerous parameters, such as rainfall, evapotranspiration, soil depth, and plant-available water. The assessment was performed using the water yield module of the InVEST model, which is mainly based on the Budyko hydrothermal coupling equilibrium assumption equation.
where Y(x) is each raster cell in the watershed x of annual water production (mm) AET(x) is the raster cellx of the actual annual evapotranspiration (mm). P(x) is the annual precipitation of the raster cell, x is the annual precipitation of the raster cell (mm), PET(x) is the annual evapotranspiration of the grid cell, and x is the potential evapotranspiration of the raster cell. Kc(x) is the crop evapotranspiration coefficient, ETo(x) is the crop evapotranspiration coefficient of the raster cell, and x is the reference crop evapotranspiration of the raster cell. AWC(x) is the effective soil water content (mm), determined by the plant-use water content (PAWC), the maximum root burial depth of the soil and the minimum value of the plant root depth. W(x) indicates the non-physical parameters of natural climate-soil properties; Z is set using the relevant literature that is available for the study area [26].

Soil retention (SDR)
The assessment of the soil retention function in the study area is based on the sediment retention module of the InVEST model, which is mainly calculated using the universal soil loss equation (USLE). The calculation equation is as follows.
where SC is the amount of soil retention t/hm 2 , which is composed of potential soil erosion Ap and actual soil erosion. Ar is determined by the difference between potential and actual soil erosion. R is the rainfall erosion factor. This paper is limited by the rainfall process data; the R value is estimated and verified according to different types of rainfall data [36]. K is the soil erodibility factor, which is calculated using the erosion-productivity evaluation model (EPIC) proposed by Williams et al. in 1990 [37]. The LS is the slope length slope factor; C is the vegetation cover and management factor; P is the soil conservation measure factor. In the model, P and C values are fixed values, which are used to reconcile the actual deviations in the calculation of soil conservation (see Table 1) and are mainly determined by referring to the USDA Handbook 703 and the related literature [38][39][40].

Net Primary Productivity (NPP)
Net primary production (NPP) refers to the net primary productivity of vegetation and is used to characterize ecosystem carbon storage services. The calculation equation is as follows. NPP where NPP(x, t) is x location, t is net primary productivity at the moment of time gc/m 2 ; APAR(x, t) is x location, t is photosynthetically active radiation at the moment of time MJ/m 2 ; ε(x, t) is x location, t is the light energy utilization at the moment of time.

Carbon Fixation Services (CF)
Ecosystem carbon stocks are composed of four main components: aboveground biomass (C above ), below-ground biomass (C below ), soil carbon (C soil ), and dead organic matter (C dead ). Aboveground biomass includes all aboveground living organic matter, such as tree trunks, branches, leaves, etc.; below-ground biomass is the living plant roots that are located underground; soil organic matter includes the organic matter fraction in the soil; dead organic matter includes dead branches and leaves, dead standing trees, etc. In addition to this, the carbon stock module of the InVEST model considers the carbon stored in the wood harvest or the associated wood products. Since this study aims to explore the status of carbon stock in the study area, the fifth carbon pool is not taken into account. The relevant equations are as follows.
C total = C above + C below + C soil + C dead (8) where C total is the total carbon stock in the study area. C above is the aboveground biomass carbon stock. C below is the belowground biomass carbon stock. C soil is the soil carbon stock. C dead is the dead organic carbon stock. The carbon density of different land-use types and different soil types in each carbon pool is mainly obtained from the relevant literature for the study area [26], where the subsurface carbon stock is converted from the aboveground carbon stock in a certain proportion.

Evaluation Index Construction of Urbanization Development Level
Urbanization is a complex system which mainly refers to the process of rural populations moving to cities, the expansion of the urban scale and the series of social changes caused by this. Based on the connotations of urban development and related research results, a relatively systematic and complete urbanization level evaluation index system was constructed. The urbanization system is divided into five dimensions, including population urbanization subsystem, economic urbanization subsystem, social urbanization subsystem and spatial urbanization subsystem, and ecological urbanization, Therefore, the evaluation of the comprehensive urbanization level (UR) needs to measure the level of the five urbanization subsystems. First, the study needs to construct an evaluation index system for urbanization development level. Demographic urbanization is the core of urbanization, and indicators reflecting the civilization of the rural population and the degree of urban population agglomeration should be considered. Economic urbanization is the driving force, and indicators reflecting regional economic development and financial income should be considered. Spatial urbanization is the guarantee, and indicators carrying human activities and guaranteeing urban development should be considered. Demographic urbanization and economic urbanization lay the foundation for urbanization and lifestyle transformation, and social urbanization lays the foundation for urbanization and lifestyle transformation. The urbanization of the population and economic urbanization lay the foundation for urbanization and lifestyle transformations; social urbanization and land-use urbanization are the expressions of society and space, respectively; ecological urbanization is an important factor in guaranteeing residents happiness, and indicators reflecting the sustainable and high-quality development of cities need to be selected. In this paper, we will follow the principles of comprehensiveness, accessibility and scientificity. We will look at the existing evaluation results of urbanization development, evaluate five different levels (population, economy, society, ecology and space), streamline and adjust their indicators, and construct an evaluation index system for urbanization level in this paper, considering development mechanisms. We will assign weights to each indicator based on the entropy value method due to their possible covariance problem, i.e., a linear relationship between the indicators. As the covariance of each index needs to be less than 7.5, the final indexes characterizing urbanization development were obtained as shown in Table 2. Exploratory Spatial Data Analysis (ESDA) is centered on measuring the degree of spatial association between things or phenomena, which compensates for the shortcomings of traditional statistical analysis, which ignores the spatially dependent properties embedded in the data [41]. ESDA analysis can be used to determine whether urbanization development has geographic clustering characteristics and spatial anomalies at ecosystem service levels, and then to explore the spatial processes in more depth. The two methods are Global Spatial Autocorrelation and Local Spatial Autocorrelation.
(1) Global spatial autocorrelation analysis is used to test whether the level of urbanization development and the level of ecosystem services are spatially correlated and spatially clustered across the study area, and is generally more commonly applied to the Global Moran's I index, with the following formula.
where S 2 = 1 n ∑ Y i − Y 2 denotes the variance in urbanization development and ecosystem service level values, n is the total number of study units. W ij is the standardized spatial adjacency weight matrix, selected based on the first-order common boundary Rook weights. When the area is adjacent to i and j, W ij = 1; otherwise, 0. Y i indicates the regional i. Y, denotes the average value of urbanization development and ecosystem service level in the study area. For I, a positive value indicates the spatial similarity of the variables, while a negative value indicates no similarity.
By the significance test of Moran's I index, if Moran's I is significant and positive at a given confidence level, it means that the level of urbanization development and the level of ecosystem services are positively correlated (higher or lower) and clustered in space, and the closer the value is to 1, the smaller the spatial difference between the level of urbanization development and the level of ecosystem services in If Moran's I is significant and negative, it indicates that the level of urbanization development and the level of ecosystem services are negatively correlated (higher or lower) between regions, and the closer the value is to −1, the greater the spatial difference between them.
(2) Local spatial autocorrelation analysis reflects the degree and significance of spatial differences between each region and its neighboring regions by measuring the significance level of spatially correlated local indicators, which is often measured by scatter plots and Local Moran's I statistics, calculated as follows: among them, w ij is the spatial weight matrix, and n is the total number of study units. By visualizing the Moran scatter plot and LISA clustering map, four major types of spatial aggregation for urbanization development level and ecosystem service level can be identified, namely, high-high aggregation, high-low aggregation, low-high aggregation, and low-low aggregation types.

Constraint Line Methods
Since the introduction of boundary lines in 1972, constraint line methods have been increasingly used in ecology to compensate for the shortcomings of traditional statistical methods in deciphering the distribution of scatter cloud data. Compared to traditional linear regression and correlation methods, constraint line methods can better characterise the limiting effects of limiting variables on response variables in complex ecosystems due to multiple factors. There are four main methods of drawing constraint lines: parametric, scatter cloud grid, quantile regression and quantile partitioning. In this paper, the quantile partitioning method was chosen to extract the constraint lines [42]. Based on the scatter plot between ecosystem services, the range of ecosystem services on the x-axis was divided into 100 equal parts to obtain 100 columns, and 99% of the quantile of each column was taken as the boundary point to obtain 100 boundary points to fit each constraint line using Origin 2021 software (www.originlab.com, accessed on 2 March 2022). The optimal constraint lines were obtained from simulations based on the shape and goodness-of-fit of the scattered point cloud to explore the constraint effects of topographic factors on ecosystem services.

Spatial Regression Model
Classical linear regression models cause the misestimation of model parameters and reduce the validity of the model when solving problems due to the presence of spatial autocorrelation; therefore, there is a need to construct regression models that are applicable to spatial data [41,43]. Anselin provided the general form of spatial regression models.
where Y is the dependent variable, X is the explanatory variable, β is the spatial regression coefficient of the explanatory variable, u is the spatially varying error term, µ is the white noise, W 1 is the spatial weight matrix reflecting the spatial decease of the dependent variable itself. W 2 is the spatial weight matrix reflecting the spatial trend of the residuals, which is usually determined based on the adjacency or distance function relationship. ρ is the coefficient of the spatial lag term, whose value ranges from 0 to 1. The closer the value is to 1, the more similar the values of the dependent variables taken in adjacent areas. λ is the coefficient of spatial error, whose value ranges from 0 to 1, and the closer the value is to 1, the more similar the values of the explanatory variables taken in adjacent regions, where W 1 can be equal to W 2 .
The following models can be derived from the general form of the spatial autoregressive model when ρ and β are not equal to 0 and λ. When these are equal to 0, the spatial lag model (SLM) is used, in which the dependent variable of the region under study is not only related to the explanatory variables in this region, but also to the dependent variables in the neighboring regions. When ρ equals 0, β and λ are not equal to 0. Therefore, a spatial error model (SEM) is used, which means that the dependent variable Y in the studied region is not only related to the explanatory variable X in the region, but is also related to the dependent variable and the explanatory variables in the neighboring regions. Spatial dependence not only implies a lack of spatial independence in observations, but also shows the structure of the data underlying this spatial correlation, that is, the strength and pattern of spatial correlation are determined by both absolute location (pattern) and relative location (distance). In this paper, we used spatial regression simulations to explore the interaction mechanisms between county-scale urbanization and ecosystem service values.

Spatial and Temporal Evolutionary Characteristics of Various Ecosystem Services in the Study Area and County
Based on the above ecosystem service evaluation methods, the spatial pattern of five types of ecosystem services at the county scale in Anhui Province were obtained and the mean values of each ecosystem service in the corresponding years were used to describe the changes in ecosystem services in the study area during the 20 years of study. The study showed that the overall performance of the five ecosystem services at the county scale in Anhui Province was "low in the north and high in the south" during the period 2000-2020 ( Figure 3). Among them, the average level of habitat support services in the study area increased and then decreased at 0.486, 0.490, and 0.485, respectively; the capacity of carbon storage services continued to decrease, with average values of 12.58 (t·ha −2 ), 12.54 (t·ha −2 ), 12.49 (t·ha −2 ). Water production service capacity increases year by year, with county average water production values of 758.07 mm, 901.87 mm and 1007.08 mm during 2000, 2010 and 2020; soil conservation function increases year by year, with county average values of 782.92 (t·ha −2 ), 782.95 (t·ha −2 ), 783.35 (t·ha −2 ). Organic matter production function first decreases and then increases, with county average values of 512.59 (gc·m −2 ), 504.98 (gc·m −2 ), 564.60 (gc·m −2 ). The overall trend of changes in ecosystem services across types showed inconsistency, mainly due to the existence of different degrees of trade-offs among ecosystem services.
(1) Analysis of synergistic relationships of ecosystem service trade-offs in the study area under different regions.
Based on the above analysis, it is worth exploring whether trade-off relationships occur among the five types of ecosystem services in the study area. Therefore, Spearman's correlation measure was used to explore the average ecosystem service correlations in different regions in Anhui Province from 2000 to 2020 at the study-area-wide scale, covering South Anhui, North Anhui, and Central Anhui (Figure 4). The study showed that the five types of ecosystem services showed strong synergistic correlations at the global scale, as well as for South Anhui and Central Anhui, while some ecosystem services showed trade-offs in North Anhui, and the overall synergistic correlations showed Anhui Province > South Anhui > Central Anhui > North Anhui. The synergistic correlation between soil conservation and carbon storage function was the best (0.75) and the synergistic relationship between soil conservation and water production was the weakest (0.58) at the whole-area scale. The synergistic relationship between habitat support services and water production services was the weakest (0.48 and 0.27) in South and Central Anhui. The water production service in this study refers to rainfall minus evapotranspiration and does not consider surface runoff, which may increase the regional water production capacity to some extent due to the increase in construction land or agricultural land. In northern Anhui, except for the strong synergy between habitat support services and carbon storage, water production services and NPP, the remaining ecosystem services showed a weak synergy and weak trade-off. The shift in the relationship between ecosystem services varied widely, e.g., soil conservation and carbon storage had the best synergy at the provincial scale, while they showed a trade-off in northern Anhui (−0.10). This is more strongly related to the different physical geographic characteristics of Anhui Province, where the northern Anhui Province is mostly plain, and its land cover/use is mostly arable land and construction land, with a smaller proportion of woodland and grassland, while the amount of woodland and grassland in the hilly areas of central Anhui and the mountainous areas of southern Anhui has a much higher proportion than that of construction land.  (2) Measurement of the effect of physical and geographical features on ecosystem service constraints.
While the relationships among the five types of ecosystem services in the study area have been discussed, whether the physical geographic characteristics of the study area are important factors affecting the trade-off relationships of ecosystem services in the study area needs to be further explored. In addition to their interactions, ecosystem service relationships are influenced by many other factors, such as meteorological factors, topographic factors, and land cover factors. The complexity of these factors also means that the relationship between ecosystem services is not simple and linear. The constraint line approach provides a new way of understanding the relationships among ecosystem services and the related ecological processes that influence them. Ecosystem patterns and processes are often influenced by a combination of factors, so a scatter plot reflecting the relationship between two ecological variables often appears as a bounded scatter cloud. Scatter clouds do not represent correlations between two variables, but instead represent constraining relationships. The constraint line approach has great potential to address the trade-offs of ecosystem services, and this section aims to explore the role of topographic factors in constraining five types of ecosystem services to better understand why relationships between ecosystems shift ( Figure 5). The results show that the constraint line of water production service according to the topography factor is "exponential"; that is, the constraint effect of topography factor on water production service continuously decreases with the rise in the topography of the whole study area (R 2 = 0.9507 The constraint line of terrain factor on soil conservation function is "convex wave type", and the constraint effect of terrain factor on soil conservation shows fluctuating characteristics (R 2 = 0.633). The topographic factor on the NPP constraint line is "exponential", and the correlation between water production service and NPP is synergistic (0.43). The two are of the same type in terms of constraint line (R 2 = 0.5908); the constraint line of topography on habitat support is "positively convex", i.e., topography effectively constrains habitat support services and proportionally increases the constraint effect, while the overall constraint line of topography on carbon storage is "convex wave" (R 2 = 0.6559). However, Figure 4 shows that the carbon storage constraint line is "open downward parabolic" at altitudes below 200 m, and the hump is low, which differs from the overall trend of soil conservation. After the elevation rises to 200 m, the constraint line fits experience a rapid rise and reach their peak at about 1000 m. The constraint fitting curve of the integrated ecosystem service function becomes "S-shaped" as the topography increases, and the constraint effect of the topography factor decreases extremely rapidly in the range of 0.6~0.7 as the integrated ecosystem service increases.

Analysis of Spatial and Temporal Changes in the Level of Urbanization Development and Subsystems in the Study Area
The population urbanization, economic urbanization, social urbanization, spatial urbanization, and ecological urbanization scores for the three period nodes of 2000, 2010, and 2020 were calculated by the weights of each index, and the spatial and temporal evolution characteristics of urbanization in the three periods were obtained from the statistical comprehensive urbanization level ( Figure 6) and divided into five levels using the natural breakpoint grading method.

Analysis of Spatial and Temporal Changes in the Level of Urbanization Development and Subsystems in the Study Area
The population urbanization, economic urbanization, social urbanization, spatial urbanization, and ecological urbanization scores for the three period nodes of 2000, 2010, and 2020 were calculated by the weights of each index, and the spatial and temporal evolution characteristics of urbanization in the three periods were obtained from the statistical comprehensive urbanization level ( Figure 6) and divided into five levels using the natural breakpoint grading method.  In the five urbanization subsystems between 2010 and 2020, the center of gravity for population urbanization gradually shifted from north to south, and spatial urbanization developed from a single point to a multi-point pattern; the pattern of economic urbanization remained the same after Hefei. The pattern of economic urbanization later maintained two major development areas around Hefei City and Wuhu City. The most dramatic change in the development level of ecological urbanization occurred in the central Anhui region, with the average level of each year rising and then falling between 2000 and 2020. In general, the comprehensive urbanization development level of counties in Anhui Province significantly increased during 2000-2020, which is in line with the normal economic and social development trends. The urbanization development in northern Anhui started to decline compared with that in central and southern Anhui, and the change in county urbanization development level levels within the region was relatively flat. The spatial pattern of development of each subsystem in Central and Southern Anhui region changed more drastically.

Measurement of the Spatial Spillover Effect of Urbanization Development on Ecosystem Services
Based on GeoDa 1.18 software (http://geodacenter.github.io/, accessed on 2 May 2021), the bivariate global autocorrelation Moran's I index between the integrated level of urbanization and the integrated level of ecosystem services in counties was measured for the years 2000, 2010, and 2020. The bivariate global autocorrelation of Moran's I index between the combined level of urbanization and the combined level of ecosystem services was −0.310, −0.312, and −0.108 in 2000, 2010, and 2020, respectively, which passed the significance test at the 0.001 level. The overall negative correlation increased and then decreased. There was a significant negative spatial correlation between urbanization level and ecosystem services during the study period, indicating that the development of urbanization leads to the deterioration of ecosystem services. The bivariate spatial autocorrelation test indicated a significant spatial dependence effect between the value of ecosystem services and urbanization, which generated a significant negative externality during the study period. This is mainly due to the rapid expansion of land for construction in the period of urbanization, which causes strong disturbance to ecosystem services. The existence of spatial spillover effects is mainly due to the flow of material, energy and information flows between neighboring regions. In addition, natural processes (e.g., atmospheric circulation, temperature, precipitation) and human activities (e.g., urbanization, economic development, population migration) also lead to the spatial spillover of ecosystem services. Therefore, ecosystem protection and management cannot be limited to a single unit, and collaborative management across regions and different levels of government departments is an effective means of ecosystem protection.
In order to further investigate the impact of urbanization development on ecosystem services, this study used a spatial regression model to investigate the impact of five types of urbanization subsystems on integrated ecosystem services (ES) to more accurately understand the impact of different urbanization types on ecosystem services. The general process of spatial regression modeling was judged with the help of Lagrange multiplier statistic (LM) in spatial econometrics. The OLS regression model was built, and the spatial autocorrelation diagnosis of the regression results was made with the help of LM-lag and LM-error (Table 3). Comparing the test statistics of LM-lag and LM-error, the spatial regression model with a better significance was selected. The results showed that both LM-lag and LM-error were significant during 2000-2020; therefore, to continue comparing the significance of Robust LM (error) and Robust LM (lag), it was found that Robust LM (error) was not significant during 2000 and 2020, and Robust LM during 2010 (error) and Robust LM (lag) were both significant, but Robust LM (lag) was more significant. Therefore, the relationship between urbanization development and ecosystem services during 2000, 2010, and 2020 was modeled using a spatial lag model (Table 4).  The study shows that the effects of population urbanization and spatial urbanization on ecosystem services during 2000-2020 are negatively correlated. After interpreting the regression coefficients measured by the SLM model, each 1% increase in population urbanization in 2000, 2010, and 2020 was shown to lead to a −1.029%, −0.341%, and −2.862% decrease in ecosystem services in that year. Each 1% increase in urbanization led to a decrease in ecosystem services of −0.140%, −0.509%, and −0.101%; social urbanization and economic urbanization led to a decrease in ecosystem services in the early stage, and each 1% increase in social urbanization will led to a −1.125% decrease in ecosystem services in that year, respectively. With social and economic development, social urbanization and economic urbanization contribute to ecosystem services. The effects of ecological urbanization on ecosystem services were all positive, with each 1% increase in ecological urbanization leading to 0.481%, 1.059%, and 0.295% increases in ecosystem services in that year. After measuring the spatial spillover effects of the five types of urbanization on ecosystem services, population concentration, economic development, and the expansion of urban construction land were found to lead to a decrease in ecosystem services. However, it is important to note that population and economic agglomeration do not necessarily lead to the degradation of ecosystem services. This is mainly due to the fact that population and economic agglomeration cannot directly interfere with ecological services, but can affect them by promoting the expansion of land for construction. The expansion of urban construction land will reduce ecological services, but inner urban waters, wetlands and green areas will mitigate the degradation of ecosystem functions to some extent; therefore, eco-urbanization promotes the improvement of ecosystem services at any stage. The spatial error term was significant during 2010 and 2020; however, in 2000 the spatial error term was not significant. This indicates that changes in ecosystem services are also influenced by other, non-urbanization factors (policies, natural environment, etc.) (Figure 7).

Transformation of the Relationship Development between Urbanization Development Patterns and Ecosystem Services in the New Era
In recent decades, the rapid development of urbanization in China has caused great disturbances to regional ecological security and seriously affected the well-being of regional people. Frequent and intense urbanization activities have led to a significant increase in land for construction, which previously had mainly shifted from arable land, forest land, and grassland, and a shift in land use driven by natural geographic features due to planning (government and residential construction needs) has begun to occur (Figure 8a). This shift involves both temporal and spatial dimensions, resulting in multiple complex adaptive changes in ecosystem services (linear, nonlinear, hysteresis, etc.). In the new era of China's national policy development shift, the development model of urbanization has shifted from the crude development of large-scale demolition and construction that occurred in the past to high-quality development. This change requires that some cities face a certain degree of actual demand for construction land, which is less than the existing construction land, resulting in the phenomenon of "hollow cities". Urbanization is the inevitable result of social development; however, in the future, there is a need to establish a long-term mechanism for the efficient, intensive and economical use of land, to fully explore the potential of urban land, to make intensive and economical use of land, to promote the more efficient allocation of resources in cities and towns, and to reduce energy consumption in the process of urbanization. In terms of urban construction patterns, there is a need to establish a compact and efficient land-use pattern to prevent excessive urban sprawl and disorderly development (Figure 8b). Under the condition of resource and energy constraints, China is facing the big challenge of high resource consumption, high energy consumption, and high carbon emissions in urbanization development. As the world's largest greenhouse-gas-emitting economy, China has made a commitment to be "peak carbon and carbon neutral". Therefore, the future development pattern of urbanization and ecosystem service enhancement needs to change from antagonistic to synergistic effects. According to existing studies, China's urbanization rate is expected to reach 75% in 2035 and 80% in 2050 according to the high-level scenario of China's urbanization development, with an additional 150-200 million people to be added in the future. This study confirms that the development of urbanization in Anhui Province has led to different degrees of decline in regional ecosystem services over the past 20 years. In terms of the development trend of the five types of urbanization subsystems, the urbanization development rate of northern Anhui may continue to lag behind that of the Hefei periphery, the Wanjiang urban belt, and southern Anhui. The large population shift to developed cities in China will cause the actual urban land demand in this area to weaken, but will increase the ecological and environmental pressure in the inflowing areas. The future urbanization development mode of the population outflow areas should focus on the transformation of the land stock to avoid the disorderly increase in urban area and waste of resources. For the county cities in southern Anhui with good economic development and an excellent ecological background, the future direction of urban construction should consider increasing the green spaces to carry out urban functions. Ecological resources such as rivers, lakes and water systems are no longer just embellish-ments for cities, but need to be transformed into important spaces that can carry industries and convert ecological values.

Conservation Strategies for Ecosystem Services in Different Regions
At present, numerous studies have shown that ecological restoration policies have positive effects on environmental quality improvements [44][45][46]. However, specific policy implementation programs need to be continuously optimized. In this paper, we investigated the nonlinear relationship between topographic factors and ecosystem services, and confirmed that the constraints of physical geography on ecosystem services are consistent with the academic understanding that natural factors, as formative factors, determine the state of the ecological background to a certain extent, and provide a quantitative basis for the implementation of future ecosystem service protection policies in the study area. The Spearman correlation was used to measure the trade-offs and synergistic relationships among ecosystem services, and it was found that the five types of ecosystem services in the study area showed synergistic relationships overall; however, some ecosystem services in northern Anhui showed trade-offs. In addition, this paper was unable to measure food production as a food supply service in the ecosystem due to the limited research data. However, related studies in this region show that there is a conflict between the ecosystem services provided by agricultural land and ecological land (woodland and grassland), which reflects the objective existence of a competitive relationship between food production and ecological conservation. Therefore, policymakers often address the possibility that ecosystem regulating services such as the conservation of regional biodiversity, air purification, and carbon sequestration and oxygen release may affect grain production. The northern Anhui region is an important supply area for grain production, so the trade-off between grain production and other ecosystem services needs to be emphasized while improving the quality of the ecological environment. The hilly area in central Anhui is the core area of development in Anhui, with high demands for urban construction and development. In the future, not only should the expansion of construction land be reasonably arranged, a certain proportion of ecological land be maintained, and low-pollution and low-consumption green industries be developed, but the improvement of internal habitat quality can be achieved by constructing ecological corridors, reducing landscape fragmentation, and increasing the connectivity of ecological land. The mountainous areas in southern Anhui Province have a high vegetation coverage, but less farmland and a smaller proportion of construction land and play an important ecological service function in Anhui Province and even in the Yangtze River Delta urban agglomeration. Therefore, attention should be paid to protecting the ecological environment in the mountainous areas of southern Anhui, and avoiding the expansion of large-scale construction land, which would block species migration corridors and destroy species diversity.

Shortcomings and Prospects
This paper measured the spatial and temporal evolution characteristics of the five types of ecosystem services in the study area at the county level during 2000, 2010, and 2020 using the InVEST model and CASA model, and explored the constraint effects of physical geographic features on the five types of ecosystems in the study area based on the constraint line method. The spatial distribution of the five types of urbanization in the county study unit was measured using the entropy method. The spatial regression model was used to measure the effects of population, economic, social, spatial, and ecological urbanization on ecosystem services, and a spatial spillover effect was found between urbanization development and ecosystem services. Additionally, various factor flows between regions were shown to affect each other in neighboring regions. Due to the difficulty of data acquisition, there are limitations in the selection of urbanization evaluation indicators in this paper, which do not fully represent the full significance of the five types of urbanization: demographic, economic, social, spatial, and ecological. Ecosystem services include provisioning services (water production), regulating services (carbon fixation, soil conservation), cultural services (spiritual, cultural and cultural benefits, etc.), and supporting services (habitat support, organic matter production), which involve many subtypes of ecosystem services. In this paper, only those ecosystem services that are closely linked to human activities in the study area were selected; this did not include plant pollination, pest and disease pests, cultural services, etc. In the future, a more comprehensive assessment of ecosystem services could be conducted using different research methods.

Conclusions
In this study, the spatial and temporal variation characteristics of the five subtypes of urbanization and ecosystem services of county units in the study area were evaluated based on the InVEST model, CASA model, constraint line method and entropy method between 2000, 2010 and 2020, and the spatial relationship between integrated ecosystem services and urbanization subtypes was further explored using spatial regression models. The main conclusions are as follows: (1) At the county scale of Anhui Province in 2000, 2010 and 2020, the overall performance of the five ecosystem services was "low in the north and high in the south". The overall trend of each type of ecosystem services was inconsistent, mainly due to the different degrees of trade-offs among ecosystem services. Using the constraint line method to measure the topography of the five ecosystem services, we found that the constraint effect of topography on the five ecosystem services is an important factor leading to the existence of different degrees of trade-offs among the ecosystem services; (2) Based on the five dimensions of population, economy, society, space, and ecology, we measured the comprehensive level of urbanization and the spatial and temporal development of subtypes in Anhui counties and found that the average level of comprehensive urbanization gradually increased (0.100, 0.187, 0.382) and the center of gravity of the highest level of county urbanization shifted from northern Anhui (2000) to central Anhui (2010) to southern Anhui (2020); (3) Based on the spatial autocorrelation model and spatial regression model, which were used to measure the spatial relationship and influence degree of integrated ecosystem services and urbanization subtypes, a significant negative spatial correlation was found between urbanization level and ecosystem services during the study period, indicating that the development of urbanization would lead to the deterioration of ecosystem services. The flow between county elements would also cause a decline in ecosystem services in the surrounding counties.