What’s the relationship between urbanization and ecosystem services from urban to rural area

: Ecosystem services (ESs) are closely bound up with human welfare. To propose regional exploitive tactics for the harmonious relationship between urbanization and ecological environment, it is essential to sort out the relationship between urbanization and ESs, and to untangle the urban-rural differences of the ecological environment effects. Regarding Beijing as an example, this study firstly explored the spatiotemporal characteristics of the ESs and the urbanization. Secondly, this study applied the latest methods based on nighttime light and K-means algorithm to identify the urban area, urban-rural fringe and rural area in Beijing. Then we explored the spatiotemporal interaction characteristics of urbanization on ESs through decoupling analysis. Finally, the correlation analysis and curve fitting were used to quantify the correlation between urbanization and ESs in urban scale and regional scale. The main conclusions included: (1) ESs declined with the increase of urbanization; (2) In cities or regions, coordinated development was the main trend of ESs and urbanization; (3) The relationship between changes in ESs and changes in urbanization has regional differences. remote sensing, GIS spatial analysis, ecological assessment models and statistical methods, “ESs - Urbanization - Urban Space Division - ESs and Urbanization" was the study clue to identify the spatiotemporal characteristics of ESs and urbanization, and to analyze the spatial interaction features and correlations between the two at urban scale and regional scale. This study had three specific purposes: (1) to use the K-means algorithm to identify Beijing’s urban area, urban-rural fringe and rural area; (2) to clarify the spatio-temporal interaction features between urbanization and ESs and grasp the trend of urbanization and ESs; (3) to explore the multi-scale relationship between urbanization and ESs. Clarifying the differences of the relationship between the same ecosystem service and urbanizations or between ESs and the same urbanization in different urban space divisions. This study can not only provide government decision-making suggestions in Beijing, but also provide the model for ecologically sustainable development in other similar regions.


Introduction
Humans derive various benefits called ESs from the ecosystem directly or indirectly (Daily, 1997), including the material resources and services that human needs, which are closely bound up with human welfare (MEA, 2003).
These benefits can be quantified by evaluating ESs (Kang et al., 2016;De Groot et al., 2012). However, the natural ecosystems have been negatively influenced by the high concentration of population, industries and the change of land use patterns, which made the ecosystem service functions decline (Peng et al., 2017).

With the process of urbanization, a large area of natural landscapes dominated by vegetation has turned into
The current research about how the urbanization affected ESs can be summarized in three aspects. The first aspect was about the variation features of the ESs with the urbanization, such as the spatio-temporal characteristics of urbanization and ESs in BTH Region (Zhou et al., 2018) as well as the trade-offs and synergy between various ESs and urbanization in Nanjing . The second aspect was the quantitative relationship between ESs and urbanization, such as the spatial correlation between ESs and urbanization (Yuan et al., 2018), correlation analysis between ESs and urbanization , or a regression model to fit the relationship between the two (Peng et al., 2016;Wan et al., 2015). The third aspect was about the influence mechanism of urbanization on ESs. It often refers to the impact of land use structure and landscape pattern with urbanization on ESs (Song and Deng, 2017;Hou et al., 2020;. Although domestic and foreign research has greatly enriched the research system related to ESs and urbanization, there are still some problems to solved. First, the relationship between ESs and urbanization is not clear enough. Most studies have shown that urbanization had a negative effect on ESs (Zhou et al., 2018), but there is still no unified conclusion about whether the relationship is linear or nonlinear (Wan et al., 2015;Song et al., 2015;Wang et al., 2019). Second, administrative division was mostly used by scholars, ignoring the urban-rural differences. Along with urbanization, cities will gradually spread to the rural-urban fringe, even rural areas, which will affect ESs in these areas. For two regions with the same degree of urbanization, the effect of natural environment on the carrying capacity for urbanization may be different (Dietz and Rosa, 1994), leading to the unbalanced spatial distribution of ESs . Thus, the impact of urbanization on ESs in urban area differed from that in rural area. However, the impact of urbanization on ESs in suburban areas have not been paid great attention to (Peng et al., 2016) or have compared the differences of urbanization on ESs between urban and rural areas.
Therefore, this study took Beijing as study area and took 2000/2005/2010/2015 as time nodes. Through the 5 remote sensing, GIS spatial analysis, ecological assessment models and statistical methods, "ESs -Urbanization -Urban Space Division -ESs and Urbanization" was the study clue to identify the spatiotemporal characteristics of ESs and urbanization, and to analyze the spatial interaction features and correlations between the two at urban scale and regional scale. This study had three specific purposes: (1) to use the K-means algorithm to identify Beijing's urban area, urban-rural fringe and rural area; (2) to clarify the spatio-temporal interaction features between urbanization and ESs and grasp the trend of urbanization and ESs; (3) to explore the multi-scale relationship between urbanization and ESs. Clarifying the differences of the relationship between the same ecosystem service and urbanizations or between ESs and the same urbanization in different urban space divisions. This study can not only provide government decision-making suggestions in Beijing, but also provide the model for ecologically sustainable development in other similar regions.

Study area
Beijing (39°28′-41°05′N and 115°25′-117°30′E) is adjacent to the Bohai Bay, lying in the northern part of North China Plain (Fig.1). In the end of 2018, the number of permanent residents reached 21.542 million in Beijing, with a per capita GDP of 21,200 dollars, which was in line with the developed countries or regions. Since the 1990s, stage of sub-urbanization was the main process in Beijing, with changes of large-scale land cover. After more than 30 years of development, Beijing's central area has surpassed the two urban greenbelts originally designated by the Beijing Municipal Government to avoid unlimited urban expansion. And now it has extended to the rural areas surrounding the central areas (Zhang et al., 2014;Zhao, 2013). At the same time, multiple satellite towns are formed around the central city and the trend is maintained. At present, it has entered the integration process of suburbanization (Peng et al., 2018), which has caused continuous reconstruction of the urban spatial pattern. Therefore, Beijing can be regarded as a typical representative of rapid urbanization areas.

Methodology Data sources
The NPP/VIIRS nighttime light data and the DMSP/OLS nighttime light data were used to differentiate urban areas, urban-rural fringes and rural areas. These two types of data are all from National Oceanic and Atmospheric Administration (NOAA) National Geophysical Data Center (https://www.ngdc.noaa.gov/eog/download.html). Then, Population grid data (https://landscan.ornl.gov/landscan-datasets), GDP grid data (http://www.resdc.cn) and urban land data (http://www.geosimulation.cn/GlobalUrbanLand.html) were used to represent urbanization of population, GDP and land. Finally, land use data (http://www.resdc.cn), annual precipitation data (http://www.resdc.cn), NPP data (http://www.resdc.cn), soil data (http://webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database/), daily weather data (http://data.cma.cn) were used to calculate the ESs. Among them, the land use data includes six primary types of cultivated land, forest land, grassland, water, construction land and unused land. The soil data includes the percentage content of clay (Clay), silt (Silt),sand (Sand),organic carbon (C) and soil reference depth (mm). The daily meteorological data includes the cumulative precipitation (mm), daily minimum 7 temperature (°C) and daily maximum temperature (°C) in Beijing.

Ecosystem service
The diverse types of ecosystems in Beijing can provide a variety of ESs for humans. At the same time, on account of natural and human disturbances, some ecological issues hinder the sustainability of Beijing (Zheng et al., 2013). Water resource of per capita in Beijing was only 4.1% of the national average in 2014. Besides, the area of soil erosion exceeded 19.5% of the total area in Beijing and the urban heat island effect was obvious.
Based on this, this study selected three typical ESs, namely water conservation, soil conservation and carbon sequestration and oxygen production, which have attracted extensive attention in the research of urban ESs (Grafius et al., 2016;Liu et al., 2017).
(1)Water Conservation Water Conservation (WC) refers to the process and ability of ecosystem to retain water in certain time and space, whose simulation was in accordance with the principle of water balance (Xu et al., 2019) and calculated using Eq. (1): where WR refers to WC (mm), P represents the average precipitation (mm), R represents the surface runoff (mm) and ET represents the surface evapotranspiration (mm). Among them, surface runoff (R) is the product of precipitation and surface runoff coefficient. Combined with literature data, the average runoff coefficients of various land use types in Beijing are obtained ( Table 1). The actual evapotranspiration (ET) is calculated according to the technical guidelines for the delineation of ecological protection red lines: where 0 represents the multi-year average potential evapotranspiration, which is calculated based on the improved Hargreaves formula (Droogers and Allen, 2002); refers to the land cover influence coefficient in the technical guidelines for the delineation of ecological protection red lines.  29 Huang, 2012;Lu, 2013;Ye, 2007Shrub 4.02 Lu, 2013 Open forest land 1. 29 Huang, 2012;Lu, 2013;Ye, 2007 Other forest land 9.14 Wang, 2011 (2)Soil Conservation Soil Conservation (SC) refers to the difference between the amount of potential soil erosion and the actual soil erosion providing there is not vegetation, whose simulation was on the basis of the revised universal soil loss equation (RUSLE) and calculated using Eq. (3): where A represents the annual average SC per unit area (t ha −1 year −1 ); R represents the rainfall erosivity (MJ mm ha −1 h −1 year −1 ), referring to the research of Rejani et al. (2016); K refers to the soil erodibility factor (t ha h ha−1 MJ−1 mm−1), calculated by the erosion-productivity impact calculator (EPIC) (Williams and Arnold, 1997); LS represents the slope length and slope factor, referring to the research of Rejani et al. (2016); C represents the coverage and management factor, which is set in accordance with the research of Hamel et al. (2015) and Mu et al. (2010); P represents the soil and water conservation measure factor, referring to research of Mu et al. (2010).
(3) Carbon Sequestration and Oxygen Production 9 Carbon Sequestration and Oxygen Production (CSOP) refers to the ability of green plants to absorb CO 2 in the air through photosynthesis and convert them into organic matter while releasing O 2 . It can not only regulate climate, but also balance the concentration of CO 2 /O 2 in the atmosphere, which has important benefits for the ecological environment. The simulation of CSOP is based on net primary productivity (NPP).
According to the photosynthesis equation, every kilogram of organic matter will absorb 1.63kg of CO 2 and release 1.2kg of O 2 . In order to ensure the comparability of CSOP in different years, the amount of CO 2 and O 2 will be converted into economic quantities. The cost of afforestation with CO 2 and O 2 is 260.9 yuan/t and 352.93 yuan/t respectively (Li & Zhou, 2016;Liu et al., 2013).

The quantity of urbanization
Urbanization is the most significant feature in the development of human society, which accompanied by economic growth (Wan et al., 2015), population increase (Zhou et al., 2018), urban land expansion (Chen, 2007) and life style changes (Bai et al., 2014). Specifically, with the development of economy, urbanization is faster and faster. Therefore, economic growth brings about increase of population and expansion of urban land, which ultimately improves people's living standards (Peng et al., 2016a). At present, four primary indicators (landscape urbanization, population urbanization, economic urbanization, and lifestyle urbanization) and several secondary indicators are accepted by many scholars to quantify the level of urbanization (Wan et al., 2015;Zhou et al., 2018). Combining previous research and taking into account the continuity, consistency and availability of data, this study selected population density (person/km 2 ), GDP density (10,000 yuan/km 2 ) and construction land ratio (%) to measure Beijing's population, economy, and land urbanization level respectively.

Urban space recognition based on night lights
According to Feng et al. (2020), based on night lights and K-means algorithm, Beijing is divided into three parts: urban zone, urban fringe zone and rural zone. In order to describe the changes of Beijing's urban space divisions from 2000 to 2015, the dynamic changes of the geographic types are summarized into three types, namely the urban infilled area (UIA), the urban expansion area (UEA) and rural area (RA). The UIA, UEA and RA accounted for 10%, 31% and 59% of Beijing, which shows that urban expansion was obvious from 2000 to 2015.

Decoupling analysis
"Decoupling" originated from physics, which means that there is no longer a relationship between two or more physical quantities . It is used to study the decoupling relationship between ecological environment and urbanization, which is of great significance for guiding the coordinated development of ecological environment and urbanization (Fang, 2018;Vehmas, 2007). This study conducts a decoupling analysis of ESs and urbanization from 2000 to 2015 to judge whether the relationship between the ESs and urbanization is coordinated. The four important parameter of the decoupling model are as follows: where ∆ represents the change of the urbanization index during the period. and −1 are the urbanization indexes at time t and t-1 respectively.
where ∆ represents the change of ESs during the period. and −1 are the ESs at time t and t-1 respectively.
where and represents the change rate of urbanization and ESs respectively during the period. Thus, the decoupling index is as follows: According to the values and relationships of the above parameters, decoupling can be divided into six categories, namely strong decoupling (SD), weak decoupling (WD), expansive negative decoupling (ER), and strong negative decoupling. (SR), weak negative decoupling (WR), declining decoupling (RD). The specific decoupling state is determined in Fig.3.

The spatio-temporal characteristics of urbanization and ESs
The spatio-temporal characteristics of urbanization From 2000 to 2015, the degree of the population urbanization, economic urbanization and land urbanization levels continuously improved. (Fig.4). Population urbanization increased the most from 2010 to 2015, rising about 38%. The growth rate of economic urbanization gradually accelerated, while the growth rate of land 12 urbanization gradually slowed down. According to spatial distribution, the three types of urbanization had similarities in each year. Population urbanization shows the characteristics of "high in the central city and low on the fringe of the city" in each year. In addition, there were differences in the distribution from year to year, which was mainly manifested in the range of high-value population density, tending to expend with time gone by. For economic urbanization and land urbanization, the spatial distribution of each year has showed the characteristics of "high in the south and low in the north, high in the middle and low on all sides". According to the change in GDP density, it can be seen that Beijing's GDP density mainly increased slightly and these regions accounted for 70% of the total area from 2000 to 2015. In addition, comparing the spatial distribution of the proportion of construction land in each year, there was no significant expansion in high-value areas of construction land, mainly because the proportion of construction land in high-value areas has further increased, especially in the areas between the fifth and sixth ring roads. According to the change, the proportion of construction land in Beijing mainly changed slightly and the area with slight decrease was wider, mainly in the ecological conservation development zone located in the north and west of Beijing. The spatio-temporal characteristics of ESs There was a significant time-space difference of WC in Beijing, while SC was similar (Fig.5). From 2000 to 2015, the WC of the ecosystem showed an upward trend, while the SC was stable. In terms of spatial distribution, WC showed the characteristics of "high in the west and north, while low in the east and south" in each year. That is, the northern and western regions of Beijing had strong WC capacity. For SC, there were two high-value regions in the north such as Yanqing, Huairou and in the southwest such as Mentougou, Fangshan.
According to the changes of WC from 2000 to 2015, there were more areas of increase than areas of decrease in Beijing. The areas grew slightly mainly, which accounted for 92.5% of the total area. There were fewer areas with slightly decrease, which scattered in northern Beijing. Areas with significant decrease and increase in WC were distributed around the periphery of the central city. Among them, the former were concentrated in the east, north and the latter were concentrated in the east, west and south.

The characteristics of spatio-temporal interaction between ESs and urbanization Analysis on the decoupling of ESs and urbanization at city scale
The decoupling results of various ESs and urbanizations from 2005 to 2015 are shown in Fig.7. It can be seen that the decoupling structure between one type of ecosystem service and different urbanization was similar.
Among them, the decoupling relationships between WC, CSOP and urbanization were all dominated by SD, while the decoupling relationships between SC and urbanizations were dominated by WD, indicating that ESs and urbanizations were mostly at the stage of high-quality coordinated development or primary coordinated development. The ratio of SD between WC and urbanizations was the largest, followed by CSOP. Therefore, the level of coordinated development between WC and urbanization was the highest. However, the decoupling relationships between the same urbanization index and different ESs were different. Such a phenomenon indicated that the influence of the same urbanization index on the different ESs had a big difference, while the effect of different urbanization indexes on the same ESs had a small difference. The decoupling relationships between WC and the three types of urbanization were similar in each urban space divisions, with SD accounting for the largest proportion. In addition, the proportion of SD in the relationship between WC and economic urbanization gradually decreased from UIA to RA, while gradually increasing in the relationship between WC and land urbanization. This showed that the capacity of WC increased with the improvement of urbanization, which made regional development more sustainable. The decoupling relationships between SC and the three types of urbanization were mainly WD in the UIA and RA, whereas WD and SD in the UEA. This showed that the capacity of SC declined with the increase of urbanization, but the growth rate of urbanization was faster than the deceleration of ESs. This further indicated that Beijing sought to develop at the expense of environment and was in a state of primary coordinated development. However, a certain proportion of SD existed in UEA. This was because these areas are close to the city center and human 16 activities were frequent, which led to great change of land use, mostly converted from cultivated land to construction land. Although the SC of construction land was weaker than that of cultivated land, land management and control measures were implemented during the process. Thus, the amount of soil erosion was less than that of cultivated land, which made SC increase in these areas. The decoupling relationships between CSOP and the three types of urbanization had obvious differences in different urban space divisions. The WD, SR and ER mainly in the UIA indicated that the urbanization and ESs developed in uncoordinated manner mainly. SD and WD were dominant in UEA, and WD, SD and ER were dominant in RA. This meant that the surrounding areas were close to the city center and undertaken the spillover of urban functions. In recent years, the development intensity had increased and the vegetation coverage on the surface had been greatly reduced, which reduced the carbon sequestration and oxygen production services in these areas. The CSOP improved in the northern and western regions of Beijing, such as Mentougou, Yanqing, Huairou, Miyun and Pinggu, caused by their good foundation of the natural environment as well as the protection measures called Guiding Opinions.

The relationship between ESs and urbanization
The spearman correlation of change in ESs and change in urbanization at the city scale and urban space divisions scale were showed in Table 2. 85% of the correlation coefficients under the urban space divisions passed the significance with a confidence of 99%, which was much higher than the 44% under city scale. In addition, by comparing the absolute value of the correlation coefficient between the two variables at the city scale and the urban space divisions scale, it was found that the absolute value of the correlation coefficients in different urban space divisions increased. It was indicated that the relationship between ecosystem service change and 17 urbanization change under the urban space divisions scale was stronger than that under the city scale, that is, there were regional differences in the relationship between ESs change and urbanization change. Therefore, it was more suitable to explore the relationship between ecosystem service change and   Table 3 showed the R 2 of curve fitting between the change in ESs and the change in urbanization in the entire Beijing area, UIA, UEA and RA. 59% of the fitting curve R 2 between the change in ESs and change in urbanization at the urban space divisions scale was larger than the fitting results at the city scale, indicating that the fitting effect of the two variables significantly improved at the urban space divisions scale. That is, the relationship between the two existed difference in different urban space divisions. Among them, the fitting curve R 2 of these four groups (SC and population urbanization, CSOP and population urbanization, SC and economic urbanization, SC and land urbanization) were all less than 0.02 and the fitting effects were not good, indicating that the relationships between these variables were not clear. In addition, WC and population urbanization, WC and economic urbanization, CSOP and land urbanization, whose fitting curve R 2 had a small gap in different urban space divisions, with the difference between the maximum and minimum of R 2 was less than 0.05 and thus the advantage of partition fitting was not outstanding. Therefore, the fitting curves of the above seven groups would not be analyzed, only discussing the fitting curves between the change in CSOP and change in economic urbanization, change in WC and change in land urbanization. The relation curves between CSOP change and economic urbanization change in UIA, UEA and RA were all quadratic functions, but the changing trends were different (Fig.9). In the UIA, there was an inverse U-shaped curve between the change in CSOP and the change in economic urbanization. When ΔGDP is 7,500, it is the inflection point of the curve. Before the inflection point, the change of CSOP increased with the increase of economic urbanization change, while the change of CSOP decreased after the inflection point. In the UEA and the RA, the relation curve between the change in CSOP and the change in urbanization was similar, all decreasing with the increase of different urbanizations change. The difference was the zero threshold for CSOP in these two divisions, that is, 58.7 million yuan/km 2 in the UEA and 44.4 million yuan/km 2 in the RA respectively. In other words, the degree of different urbanization in UEA would increase more when the CSOP didn't decrease. In addition, the change in CSOP caused by change in economic urbanization were not the same in different divisions. The change in economic urbanization within the UIA corresponded to smaller changes in CSOP. This is because the degree of urbanization in UIA was high and the ecosystem service foundation was weak. Therefore, more attention was paid to the protection of the ecosystem during economic development, such as the construction of road green belts, parks, etc., which weakened the negative effects of urbanization on the ecosystem to a certain extent. However, there was more urgent needs for urban expansion in the UEA and RA, and thus the lack of reasonable planning for ecological protection in the process of urbanization led to significant damage to ESs.
There is a large gap in the fitting effect of the relation curve between change in WC and change in land urbanization in the three urban space divisions. Among them, the fitting R 2 in UIA and RA were less than 0.01, indicating that the relationship between the variables is not clear. In the UEA, the fitting curve R 2 of the change in WC and change in land urbanization was 0.0787, which was much higher than the other two divisions, suggesting that the curve fitting between the change in WC and the change in land urbanization was disturbed by the scattered distribution in UIA and RA. In the UEA, there was an inverse U-shaped relation curve between the WC change and the land urbanization change. When Δclp is 0.01, it is the inflection point of the curve. Before the inflection point, the WC change increased with the urbanization change increased, while the WC change 20 decreases after the inflection point. It was worth noting that when the increase in the proportion of construction land was more than 0.26, the change in WC became negative. compared with that in 2000, the WC improved in most streets. These streets were distributed evenly in various urban space divisions, so the decoupling relationships between WC and urbanization appears insensitive to the urban space divisions. Second, although human activity interference and land use changes were serious in UEA than that in UIA and RA, land management and control measures were implemented in the process. Thus, the amount of soil erosion was less than the original, and SC increased instead. This explained that the proportion of streets with strong decoupling relationship between SC and urbanization in the urbanization expansion zone was much larger than that of the other two divisions. It also explained the proportion of decoupling types in the UEA and the non-UEA. Finally, due to the vegetation type and cover rate, the NPP production was obviously different in different urban space divisions. Therefore, there were also obvious differences of the decoupling relationships between CSOP and urbanization in different divisions. To sum up, the sensitivity of the decoupling relationship to divisions lies in the sensitivity of each ecosystem service to land cover. CSOP completely depends on the vegetation type, whereas SC is the second most affected by the vegetation type, followed by WC.
Many scholars have paid attention to the ESs in different urban space divisions. For example, the research results of Xu et al. (2020) and Hou et al. (2020) were similar to our study, which showed that ESs gradually decreased from urban centers to rural areas. However, few scholars have studied the relationship between ESs and urbanization in different urban space divisions, and most of them have analyzed the relationship between ESs and urbanization on the city scale. CSOP and economic urbanization were in an inverse U-shaped curve in UIA, while there was a negative correlation in UEA and RA. However, the growth rate of CSOP per unit of economic urbanization in UEA was larger than that in RA. The reason may be the level of urbanization was high in the UIA and the ecosystem service foundation was weak. Therefore, more attention is paid to the protection of the ecosystem. More protection measures will be taken during development, such as the construction of road 22 green belts, parks, etc. These measures have weakened the negative effects of urbanization on the ecosystem to a certain extent. However, needs for urban expansion in UEA and RA was more urgent. The lack of reasonable planning for ecological protection in the process of urbanization has caused significant damage to ESs. The relationship between WC and land urbanization changes is most obvious in urban expansion areas, showing an inverse U-shaped relationship. The possible reason was that with the help of the relevant land protection policy, the increase in the proportion of construction land would not affect the WC of Beijing within a certain threshold, but illegal land occupation would occur after a certain threshold, leading to the decline of WC.
The innovation of this study is to explore the relationships between change in ESs and change in urbanization based on the three types of urban space divisions instead of the urban-rural dual structure or city scale of the previous study. Through decoupling analysis and curve fitting, the relationship between ESs and urbanization can be clear more. However, there are some shortages as follow in this study. First, the types of ESs are limited and more types should be included in the study. Based on the actual situation in Beijing, this study focused on three typical ESs, namely WC, SC, and CSOP, but lacks important ESs such as food supply and cultural services. As a megacity and the capital of five dynasties, the above-mentioned services are all important components of Beijing's urban ecosystem and need to be considered in the next step of research. Second, the research results will be affected by the granularity of the analysis unit and grid-based correlation analysis can be carried out to explore the differences in the relationship between ESs and urbanization based on different granularities. Finally, the mechanism of relationships between urbanization and ESs need to be studied. This study initially explored the relationship between ESs and urbanization at the city and divisions scales, but the internal mechanism of the two needs to be further studied. In the future, combined with the distribution characteristics of the data, spatial panel models such as spatial lag models and spatial error models can be used to explore the driving factors of urbanization affecting ESs, and to assess the direct and indirect effects of independent variables on dependent variables.

Conclusions
(1) ESs decline with the increase of urbanization. The level of urbanization continued to improve from 2000 to 2015, and the ESs, in addition to the relatively stable SC, WC and CSOP increases. In terms of spatial 23 distribution, the urbanization level and ecosystem service capacity show opposite trends from the central city to rural area, that is, the urbanization level gradually decreases and the ecosystem service capacity gradually increases.
(2) Coordinated development is the time-space interaction characteristic of ESs and urbanization. In Beijing, WC, CSOP, and urbanization are in a state of high-quality coordinated development in Beijing, while SC and urbanization are mainly in a primary coordinated development state. Among them, WC and urbanization in the urban infill areas and rural areas are high-quality coordinated development, SC, CSOP and urbanization are primary coordinated developments; urban expansion areas serve various ecosystems and urbanization are high-quality Coordinated development.
(3) The relationship between changes in ESs and changes in urbanization has regional differences. CSOP and economic urbanization are in an inverse U-shaped relationship in UIA, while there is a negative correlation in urban expansion areas and rural reserved areas. However, the growth rate of CSOP per unit of economic urbanization in urban extension area is higher than that in rural areas. The relationship between WC and land urbanization changes is most obvious in urban expansion area, showing an inverse U-shaped relationship.

Suggestions
(1) Along with urbanized level enhancement, the increase of the population density in urban expansion area is faster than that in urban infilled area and rural area. This is due to the fact that housing prices in urban expansion areas are lower than that in core areas. Besides, the infrastructure in urban expansion areas is also better than that in rural areas, which is attractive for many people. But it will affect negatively ESs if the population exceeds a certain threshold. Therefore, urban New constructions can be built to evacuate the population and industries in the central area, and be equipped with necessary facilities. To avoid urbanization at the expense of the environment, the environmental protection line should be supported to protect natural resources before development. In this way can the urban economy speed the development, but the environment is also better and better. Besides, urban with infill expansion instead of blind expansion is preferable.
(2) The level of urbanization in urban infilled area is high, but the ecosystem service foundation is weak.
Therefore, more attention should be paid to the protection of the ecosystem and more protection measures should be taken, such as the construction of road green belts and parks. These measures have weakened the negative 24 effects of urbanization on the ecosystem to certain extent. Because need for urban expansion in urban expansion area and rural area is urgent, there is a lack of reasonable planning for ecological protection in the process of urbanization, which leads to significant damage to ESs. Therefore, protection measures should be into consideration in these areas to avoid protecting after destroying.
(3) In order to achieve a coordinated development between CSOP and urbanization, the increase of GDP density in urban infilled area needs to be within 7,500 so that CSOP increases with urbanization rising. However, the rapid speed of economic urbanization needs to be moderately controlled in urban expansion area and rural area and the reserve of CSOP should be protected. WC in urban expansion area is influenced by the glass in construction land. The increase of the construction land density in urban expansion area needs to be within the range of 0.01 so that the urbanization can be in harmony with carbon and release oxygen production.

Declarations
Funding: This study was financially supported by the National Key Research and Development Program of China (2019YFB2102000).

Competing interests:
The authors declare no conflicts of interest.

Availability of data and materials:
The datasets used or analyzed during the current study are available from the corresponding author on reasonable request.