Spatially Heterogeneity Response of Critical Ecosystem Service Capacity to Address Regional Development Risks to Rapid Urbanization: The Case of Beijing-Tianjin-Hebei Urban Agglomeration in China

Abstract

Urban agglomerations have become the new spatial unit of global economic competition. The intense socioeconomic activities attributed to the development of urban agglomerations are bound to cause damage to the ecosystem services of these urban agglomerations. This study adopts the Beijing-Tianjin-Hebei urban agglomeration in China as the research object, analyzes the spatiotemporal evolution of its critical ecosystem service capacity to address regional ++-development risks from 2000–2018, and employs the Moran’s I and geographically weighted regression model to explore the spatial correlation and spatial heterogeneity in the responses of urbanization and ecosystem services. The study indicates that (1) from 2000–2018, the ecosystem services of the Beijing-Tianjin-Hebei urban agglomeration exhibit an increase and then a decline, reaching the highest index in 2015; (2) the ecosystem services reveal obvious spatial heterogeneity with the Yan and Taihang Mountains region as the boundary; (3) built-up area ratio, GDP density, and population density exhibit highly obvious negative correlation driving characteristics on ecosystem services; and (4) the construction land ratio exerts a notable impact on areas with a high ecosystem services, while the spatial response of the effect magnitude of the population and GDP densities is largely influenced by intensive, high-pollution and energy-consuming industries. This article also proposes strategies for the optimization of ecological resources and spatial control, which are dedicated to mitigating the negative impacts of rapid urbanization processes on ecosystem services.

1. Introduction

Since the beginning of the 21st century, urban agglomerations have increasingly become the new spatial unit of global economic competition due to rapid urbanization and expansive population growth, thereby increasing globalization and intensifying competition on a global scale [1,2,3]. According to United Nations projections, by 2050, 75% of the world’s total population will live in cities [4], and the population of the 40 largest urban agglomerations or metropolises will account for 18% of the global total population. These cities or metropolises host approximately 66% of the global economic activity and 85% of the total technological innovations [5]. As areas of rapidly growing economic, urbanization and industrialization, the development of urban agglomerations could inevitably lead to changes in land use patterns and their extent, which has a considerable impact on natural resources, ecosystem functions and biodiversity, leading to ecological crises such as freshwater scarcity, air pollution and regional climate change [6,7,8], posing a serious threat to regional ecosystems [9,10,11] and affecting the ability of ecosystems to provide services [12,13]. Ecosystem degradation caused by rapid development of urban agglomerations has become a global problem with drastic impacts on human well-being and economic development [14,15]. This problem could become increasingly acute in the coming decades [16]. Therefore, clarification of the impact of urbanization on the ecosystem services (hereafter referred to as ESs) of urban agglomerations is of critical importance to achieve sustainable regional development.

Currently, urban agglomerations are facing the challenge of rebalancing ecosystem conservation and urbanization to maximize the benefits for the region. ESs are tangible and intangible benefits provided by ecosystems for human social development, which can provide many services such as carbon fixation, fresh water, and air improvement. Some scholars have begun to conduct research on integrated regional development from an ecological perspective. Exploring the relationship between the impact of urbanization on ESs has become mainstream in providing regional ecosystem management and land use planning decisions. Many studies have shown that urbanization has a significant impact on the ESs, and it has become increasingly difficult to reconcile the growing demand for urbanization with ecological sustainability. On the one hand, urban expansion converts ecologically capable woodlands and grasslands into land for construction, destabilizing ecosystems [17]. On the other hand, the surrounding areas provide a large amount of natural resources for urban development, and the large amount of pollution emissions caused by urban development also poses a threat to regional ecological security [18,19].

Faced with the increasing problem of regional ecological risks due to urbanization, scholars have conducted many studies on how ESs will be affected by urbanization, although most studies urbanization on ESs have focused on the overall or single type of service impacts of urbanization on the value of Ess [20,21,22]. However, a minority of scholars have begun to focus on the impacts of urbanization on critical ESs that maintain regional ecological security. It was found that rapid urban expansion in China often triggers changes in the regional natural environment [23], causing a series of ecological problems such as regional soil erosion, low-lying water bodies, reduced carbon storage, soil erosion and loss of biodiversity, weakening ESs and threatening regional sustainable development [24,25,26], with biodiversity conservation, outdoor recreation, water production, carbon storage, heat island mitigation ESs such as biodiversity conservation, outdoor recreation, water production, carbon storage, heat island mitigation and soil erosion are the most sensitive to urbanization [27,28,29,30]. Overall, previous studies have focused more on the impact of urbanization on which ESs that maintain regional ecological security, but there is still a lack of research on the spatial heterogeneity of the impact of urbanization on critical ESs to guide the synergistic development of regional urbanization development and ecosystem conservation [31,32].

China is one of the most rapidly urbanizing regions in the world [33,34]. The Beijing-Tianjin-Hebei urban agglomeration (hereafter referred to as the BTH) represents a major national strategy in China, and the development goal of forming a world-class large urban agglomeration with Beijing and Tianjin as the core cities has been established. The BTH is the most densely populated region in northern China and one of the regions with rapid urbanization in China. This paper selects the core region of the BTH as the research object. In order to avoid ignoring the rich information inherent in ESs at spatial and temporal scales in the calculation of ESs using Land use and land cover, we use the ES index to evaluate Ess, and constructing Moran’s I and geographically weighted regression model (hereafter referred to as the GWR model) to explore the spatial heterogeneity of the impact of urbanization on critical ESs. This objectives of this study are: (1) construct an indicator system for critical ESs in the BTH; (2) quantitatively analyze the spatial and temporal correlations between urbanization processes and critical ESs; and (3) determine the spatial heterogeneity of ESs responses to urbanization processes. In addition, we attempt to answering the following questions: (1) what spatial and temporal dynamic evolutionary characteristics do critical ESs in the BTH exhibit? (2) How do the urbanization processes of each type in the BTH change the spatial and temporal dynamics of critical ESs in the region? (3) How to synergize regional urbanization development and ecosystem conservation in the future?

2. Materials and Methods

2.1. Study Area and Study Data

In this study, eight cities, namely, Beijing, Tianjin, Chengde, Tangshan, Zhangjiakou, Baoding, Cangzhou, and Langfang, were delineated as the BTH city cluster study area (Figure 1), located at 113°40′–119°19′ E longitude and 37°29′–42°37′ N latitude, with a total area of 160,000 km². The 2018 resident population of the region amounted to 74.868 million people, and the GDP in 2018 reached up to 680.8146 billion yuan.

The study area is located north of the lower reaches of the Yellow River, with the Yanshan Mountains to the north, the North China Plain to the south, the Taihang Mountains to the west, and Bohai Bay to the east, with a typical temperate monsoon climate. The topography gradually slopes from northwest to southeast, precipitation decreases from southeast to northwest, and there occur various types of landforms and ecosystems, which are divided into five major regions from northwest to southeast, namely, the Bashang Plateau Region, Yanshan and Taihang Mountains, Beijing-Tianjin-Bao Center Area, Coastal Sea Region, and Plain Region. According to 2018 statistics, the average NDVI in this region was 0.753, in which the arable land area accounted for 45.87%, the forestland area accounted for 22.36%, the grassland area accounted for 16.80%, and the water area accounted for 3.55% of the total land area. The construction land area reached 16,876.51 km², accounting for 10.50% of the total area, and the urbanization of BTH increased by 12.1% in a decade, which typifies urbanization development in northern China. The urbanization rate reached as high as 70%, and the study area is the most rapidly developing urbanization region in northern China, which occupies an important strategic position in the economic and social development of China and exhibits a high potential to develop into a world-class city cluster.

The research data include seven categories: basic geographic information data, land use data, remote sensing (NNP and NDVI) data, digital elevation model (DEM) data, meteorological data, soil property data, and socioeconomic statistics. Basic geographic information data and administrative divisions were obtained from the Data Center for Resource and Environment Science of the Chinese Academy of Sciences (www.resdc.cn (accessed on 6 August 2021)). Land use data were also obtained from the Data Center for Resource and Environment Science of the Chinese Academy of Sciences (www.resdc.cn) and were divided into 6 categories: arable land, forestland, grassland, water body, construction land, and unused land. Remote sensing data and elevation data were obtained from the National Geographic Information Spatial Data Cloud platform (http://www.gscloud.cn/ (accessed on 12 February 2021)). Meteorological data such as temperature, rainfall and soil property data were acquired from the Data Center for Resource and Environmental Sciences of the Chinese Academy of Sciences (www.resdc.cn). The resolution of the above data is 1 km × 1 km. 2000/2005/2015 GDP and population data were retrieved from the Data Center for Resource and Environmental Sciences, Chinese Academy of Sciences (www.resdc.cn), while 2018 GDP and population data were obtained from the 2019 Statistical Yearbook of Beijing, Tianjin, and Hebei cities and counties.

2.2. Research Framework

As shown in Figure 2, first, in response to the current risks faced across the BTH, such as increased carbon emissions, high temperatures, heavy rainfall and flooding, air pollution, and soil erosion [35,36,37,38], we selected representative evaluation factors of ESs, adopted a dimensionality reduction approach to characterize ESs, and analyzed the spatiotemporal variation in ESs. Consequently, the global bivariate Moran’s I was adopted to evaluate the correlation between urbanization and ESs. The local bivariate Moran’s I model was applied to analyze the spatial aggregation phenomenon of urbanization and ESs in the surrounding districts and counties. The GWR model was applied to analyzing the spatial heterogeneity of ecosystem service responses to urbanization. Finally, through analysis of the obtained results, we proposed ecological optimization and ecosystem management measures and policies for the BTH.

2.3. Ecosystem Services and Urbanization Measurements

2.3.1. Construction and Measurement of the Ecological Services Index System

Based on previous studies of ecological services in relevant regions, three criteria were adopted to determine ESs indicators of the study area [39,40,41,42]: (1) each selected ecological services evaluation factor was consistent with the type of ESs representing the Millennium Ecosystem Assessment System; (2) the selected indicators were consistent with the demand for the various types of ESs in the BTH; and (3) relevant data were accessible. According to the current problems of soil erosion, severe urban heat island effect, weak carbon sequestration capacity, vegetation degradation, poor biodiversity, and uneven distribution of recreational services in the BTH [43,44,45,46,47,48], an index system for ESs evaluation in the BTH study area was constructed, thereby identifying one provisioning services (net primary productivity (NPP)), four regulating services (soil and water conservation, carbon and oxygen sequestration, water conservation, and heat island mitigation), one supporting services (biodiversity conservation), and one cultural services (recreation suitability) (

Table 1. Ecosystem services indicators.

Ecosystem Services Classification	Ecosystem Services Types	The Significance of Metrics	Formula		The Descriptions of Metrics
Provisioning service	Net primary productivity (NPP)	It is the fixed energy or organic matter produced per unit area and unit time remaining after green plant respiration	NPP = GPP − Ra	(1)	where Ra denotes the consumption of autotrophic respiration and GPP denotes the total primary productivity.
Regulating services	Soil and water conservation	It is the ecosystem function preventing soil erosion caused by water erosion through corresponding structures and processes.	Spro = NPPmean × (1 − K) × (1 − Fslo)	(2)	where Spro is the soil and water conservation serviceability index, Fslo is the slope factor, and K is the soil erodibility factor.
	Carbon sequestration and oxygen release	It represents material exchange processes involving oxygen and carbon dioxide through photosynthesis by green plants.	Gv = 1.63Rc × A × NPP; G0 = 1.19A × NPP	(3)	where Gv is the annual carbon sequestration by vegetation, G0 is the annual oxygen release, Rc is the carbon content in CO2.
	Water conservation	It suggests that an ecosystem interacts with water through its unique structure to intercept, infiltrate, store, regulate water flow, the water cycle through evapotranspiration.	WR = NPPmean × Fsic × Fpre × (1 − Fslo)	(4)	where WR is the ecological spatial water holding capacity services index, Fsic is the soil infiltration factor, Fpre is the average precipitation factor, and Fslo is the slope factor.
	Heat island mitigation	It is a phenomenon of urban microclimate change in the local temperature, humidity, and air contamination across the urban surface caused by anthropogenic changes.	Ts = Ti + A(Ti − Tj) × B	(5)	where Ts is the surface temperature, Ti and Tj are the brightness temperatures in thermal channels i and j, respectively, and A and B are coefficients determined by factors.
Supporting service	Biodiversity maintenance	It which aims to preserve the roles of various species and the corresponding genetic diversity.	Sblo = NPPmean × Fpre × Ftem × (1 − Falt)	(6)	Where Sblo is the biodiversity maintenance services capacity index, Fpre is the mean rainfall, Ftem is the mean air temperature, and Falt is the elevation factor.
Cultural service	Recreation suitability	It was assessed through including evaluation of recreation resources, recreation facilities, and recreation areas three aspects to determine whether the use of recreation resources is reasonable.	$S={\displaystyle \sum }_i^n{W}_i\times X_i$	(7)	where S is the comprehensive evaluation index of recreation use, Wi is the rank value of recreation factor i, Xi is the weight value of the different recreation factors, and n is the number of evaluation factors.

).

NDVI varies with the greening rate, which is related to a variety of important ecosystem functions, such as plant productivity, biomass, biodiversity, and water containment and carbon storage, and NDVI can reflect the overall changes in ESs more accurately [49,50]. Therefore, in this paper, first, by comparing mean NDVI from 1998–2018, four representative time points—2000, 2005, 2015, and 2018—were selected based on mutation points and combined with the development history of the BTH [51].

2.3.2. Comprehensive ES Index

We use the ES index to assess regional ESs. The following steps are involved in the calculation process [52]. (1) This calculation starts with a maximum-minimum normalization of the 7 service indicators mentioned above. (2) Accumulation is conducted according to the 4 major ecosystem service types and normalization is done here to accumulate the composite ecosystem service index.

$${\mathrm{x}}_{\mathrm{ij}}=\frac{{\mathrm{x}}_{\mathrm{ij}}-{\mathrm{x}}_{\min }}{{\mathrm{x}}_{\max }-{\mathrm{x}}_{\min }}$$

(8)

where y_ij and x_ij are the standardized and raw data of ecosystem service factor i; x_max and x_min are the maximum and minimum values of raw data of state index j.

$$\mathrm{CIES}={\displaystyle \sum }_{i=1}^n{ES}_i$$

(9)

CIES is the composite index of ecosystem services; ES is the standardized result of the i-th category of ecosystem services; n is 4.

2.3.3. Selection of Factors to Reflect the Level of Urbanization

The impact of urbanization on ESs encompasses a dynamic process involving social, economic and geographical factors and exhibits complex mechanisms of action [53]. Exploring the impact of urbanization on ESs could help us better manage the relationship between humans and nature. Therefore, in this study, three indicators widely applied to measure the level of urbanization, namely, population growth, gross domestic product, and built-up area ratio, were selected to characterize the level of urbanization in the BTH and explore their relationship with Ess [54,55,56]. Among these indicators, the population density (POPD, person/km²) was considered to reflect population urbanization, the GDP density (GDPD, CNY/km²) was employed to describe economic urbanization, and the built-up area ratio (CLP, %) was applied to indicate spatial urbanization [57,58].

2.4. Analysis Method

2.4.1. Spatial Correlation Test

Exploring whether an attribute value exhibits characteristics such as spatial aggregation of distribution and spatial anomalies is the core of exploratory spatial data analysis methods. Global Moran’s I and local Moran’s I (LISA) are spatial analysis models that measure and test spatial convergence or heterogeneous association patterns to reveal the spatial distribution characteristics of objects [59,60]. Global bivariate Moran’s I analyzes the spatial correlation between urbanization and ESs across the entire study area. The global Moran’s I value is generally between −1 and 1. If it is greater than 0, it means that the space is positively autocorrelated; if it is less than 0, it means that the spatial entities are discrete; if its value tends to 0, it means that the space obeys random distribution, and the results were subjected to Monte Carlo tests and 999 randomized simulations (p ≤ 0.05, |Z| > 2.58). The formulas used to calculate Global bivariate Moran’s I are as follows:

$$I=\frac{n}{{{\displaystyle \sum }}_{i=1}^n{{\displaystyle \sum }}_{j=1}^n{W}_{ij}}\ast \frac{{{\displaystyle \sum }}_{i=1}^n{{\displaystyle \sum }}_{j=1}^n{W}_{ij}(x_i-\overline{x})(x_j-\overline{x})}{{{\displaystyle \sum }}_{i=1}^n{(x_i-\overline{x})}^2}$$

(10)

where I = the global Moran’s I index; n = the number of spatial cells; x_i, x_j = the observed values of cell i and cell j, respectively; W_ij is the spatial weight adjacency matrix of spatial cell i and j (ij = 1, 2, 3, … n).

Bivariate LISA produces output clustering maps that illustrate the relationship between the ESs value of a particular location and the mean urbanization level of neighboring locations at a certain significance level, providing a more intuitive understanding of the local spatial correlation [61]. The final generated LISA clustering maps revealed four types of local spatial autocorrelation: the high–high (H-H) high-ES surrounded by high-urbanization; the high–low (H-L) high-urbanization surrounded by low-ES; the low–high (L-H) low urbanization values surrounded by high-ES; and the low–low (L-L) low urbanization values surrounded by low-ES, which is used to illustrate the spatial relationship between the ESs of a particular site and the level of urbanization of neighboring site. The formulas used to calculate Moran’s I are as follows:

$$I_{kl}^i=z_k^i{\displaystyle \sum }_{j=0}^n{w}_{ij}{z}_l^j$$

(11)

where $z_k^i=\frac{x_k^i-{\overline{x}}_k}{{\lambda }_k}$, $z_l^i=\frac{x_l^i-{\overline{x}}_l}{{\lambda }_l}$; w_ij is the spatial connection matrix between spatial units i, j; $x_k^i$ is the value of attribute k of spatial unit i; $x_l^i$ is the value of attribute l of spatial unit j; ${\overline{x}}_k$, ${\overline{x}}_l$ are the mean values of attributes k, l; λ_k, λ_l are the variances of attributes k, l.

2.4.2. Geographically Weighted Regression (GWR)

Studying the heterogeneity of the impact of drivers on the study object in different regional spaces and their complexity helps to deeply assess the relationship between the two in local areas [62,63]. THE GWR model is a spatial analysis technique capable of reflecting the non-stationarity of parameters in different spaces, compared to the traditional models, which consider the local impact of spatial objects and reflect the non-stationarity of parameters in different spaces, showing the neglected local characteristics [64,65,66,67]. In this paper, GWR 4.0 software is used to analyze the degree of local response of ESs to different types of urbanization at the county scale of the BTH. Meanwhile, the GWR model requires a geographic variability test for the variable coefficient independent variables, which is useful for verifying whether the theoretical assumptions are reasonable. When the DIFF of Criterion is positive, it means that the hypothetical premise of the variable coefficient term of the model is rejected and the impact of the variables is consistent at the global scale, and vice versa, it means that the variable coefficient hypothesis test is passed and the impact of the variables is considered to have a regional variability effect. The formulas used to calculate GWR model are as follows:

$$\begin{gathered} y_i={\beta }_0\left(U_i,V_i\right)+{{\displaystyle \sum }}_{k=1}^n{\beta }_k\left(U_i,V_i\right)xk\left(U_i,V_i\right)+{\epsilon }_i, \\ \beta k\left(U_i,V_i\right)={(X^{T}W\left(U_i,V_i\right)X)}^{-1}{X}^{T}W\left(U_i,V_i\right)y \\ {W}_{ij}=\exp \left(-\frac{d_{ij}^2}{h^2}\right) \end{gathered}$$

(12)

where β₀(U_i, V_i) is the geographically weighted regression intercept of the (U_i, V_i) spatial location, β_k(U_i, V_i) is the weighted regression coefficient of the kth independent variable (driver) at the (U_i, V_i) spatial location, xk(U_i, V_i) is the value of the kth independent variable (driver) taken at the (U_i, V_i) spatial location, ε_i is the residual of the algorithm, X^T is the independent variable (driver factor) transpose, W(U_i, V_i) is the distance weight matrix, and h is the bandwidth of the AIC criterion, which is the distance between spatial position i and spatial position j.

3. Results

3.1. Spatiotemporal Evolution of ESs

The evolution of the ES index in the study area is shown in Figure 3. First, the provisioning service exhibited a steady increase, from 0.241 to 0.310, during the study period. The largest increase of 0.06 occurred from 2000 to 2005. The growth trend substantially declined after 2005, reaching a maximum index of 0.316 in 2015. Second, the regulating services reached the maximum index of 0.333 in 2015, before decreasing to 0.212 by 2018. Third, the supporting service fluctuated the most among the four ESs during the study period, increased from 0.038 in 2000 to a maximum index of 0.057 in 2005 and then decreased to a minimum index of 0.020 in 2015 and increased to 0.048 in 2018. Fourth, the cultural service varied slightly throughout the study period, reaching a maximum index of 0.579 in 2005.

The spatial characteristics of the evolution of the ESs are specified as follows: (1) the provisioning service indicated differential development characteristics during the study period, with the Yan and Taihang Mountains Region as its boundaries. The Beijing-Tianjin-Bao Central Area, Yan and Taihang Mountains Region, and the Plain Region generally exhibited an increasing trend. In contrast, the Bashang Plateau Region exhibited a trend of increasing and then decreasing during the study period. (2) In terms of the regulating services, the Bashang Plateau Region indicated a trend of annual weakening during the study period. In the Yan and Taihang Mountains Region, the ES in the northern Yan Mountains area gradually increased, but the western Taihang Mountains exhibited the increase-decrease-increase trend. The ES the Beijing-Tianjin-Bao Central Area increased annually. The Plain Region were enhanced before 2015 and then exhibited a decreasing trend. (3) The supporting service revealed a low index and did not substantially change during the study period, with the exception of the Yan and Taihang Mountains. The Yan and Taihang Mountain area exhibited fluctuations involving strengthening, weakening, and strengthening. (4) The cultural service was enhanced before 2005 and then weakened in the Yanshan and Taihang Mountain areas, the Beijing-Tianjin-Bao Central Area, and the Plain Region. The Coastal Sea Region exhibited the reverse trend. In contrast, the Bashang Plateau Region did not change much during the study period (Figure 4).

The changes in the Comprehensive ES index over the study period are shown in Figure 5. The value of ES index in the study area showed an increasing and then decreasing trend in time variation, with the highest index of 1.890 in 2015. The spatial characteristics of the ES index evolution in Beijing-Tianjin-Bao central area, plain area, Bashang Plateau area and Taihang mountain area follow the trend of “Increase-Decrease”. However, the ES index in the Yan Mountains Region exhibited a “Increase-Decrease-Increase” trend.

3.2. Characteristics of the Spatiotemporal Evolution of Urbanization

The results are shown that the three high-urbanization factor types were concentrated in the plains, while in the Yan and Taihang Mountains and the Bashang Plateau Region, there occurred significantly less urbanization. The change of urbanization level from 2000 to 2018 in Figure 6, which including POP, GDP, and CLP per unit area of BTH. In terms of the spatial expansion of urbanization, the medium–high and high values of the construction land share were mainly concentrated in the central urban areas of Beijing and Tianjin throughout the study period. The high value areas in the central urban areas of Beijing and Tianjin did not change significantly during the study period, and the main growth areas of urbanization were concentrated in Beijing, around the central urban areas of Tianjin and in the central urban areas of the remaining cities. Under the influence of urbanization agglomeration effect, the GDP density values gradually converged to the central urban areas of Beijing and Tianjin during the study period from 2000 to 2018. In addition, the GDP density showed an overall increasing trend in the study area during 2015–2018. In terms of the population density, the more densely populated areas were clustered in the plain regions, and the high value areas of population urbanization are concentrated in the central urban areas of Beijing and Tianjin. While between 2005 and 2018, the population gradually clustered in the regional space of the central urban areas of Beijing, Tianjin, and prefecture-level cities.

3.3. Effect of Urbanization on ESs

3.3.1. Global Moran’s I Analysis

This paper applied the GeoDa spatial analysis tool to create a spatial weight matrix for each urbanization indicator and E binary spatial autocorrelation analysis by the district and county, and the correlation between urbanization and ESs was evaluated through Monte Carlo simulations (ranking: 999; p value: 0.01; absolute Z value > 2.58). Global bivariate Moran’s I results indicated a significant negative spatial correlation between urbanization and ESs (all Moran’s I values < 0; p value = 0.01) (

Table 2. Moran’s I between the ESs and the GDP density, POP density and CLP.

	Variables	2000	2005	2015	2018
Moran’s I	CLP	−0.537	−0.557	−0.568	−0.576
	GDPD	−0.519	−0.479	−0.291	−0.407
	POPD	−0.488	−0.384	−0.549	−0.520
p value	CLP	0.001	0.001	0.001	0.001
	GDPD	0.001	0.001	0.001	0.001
	POPD	0.001	0.001	0.001	0.001
Z value	CLP	−10.985	−11.414	−11.483	−11.667
	GDPD	−10.730	−10.23	−7.019	−9.070
	POPD	−10.210	−8.606	−11.449	−10.953

). After 2000, the strength of the relationship between the effects of urbanization on the ESs followed, in descending order, land urbanization (CLP, Moran’s I = −0.529), economic urbanization (GDPD, Moran’s I = −0.529) and population urbanization (POPD, Moran’s I = −0.529). In 2005, there was a significant decrease in the intensity of the impact of economic and population urbanization on the ESs (GDPD, Moran’s I: −0.479; POP, Moran’s I: −0.052), except for an increase in the intensity of the impact of land urbanization on the ESs (CLPD, Moran’s I: −0.556). In 2015, this period exhibited a significant decrease in the intensity of the impact of economic urbanization on the ESs (GDPD, Moran’s I = −0.291) and a significant increase in the impact of population urbanization on the ESs (POPD, Moran’s I = −0.549), which undoubtedly led to a change in the ranking of the intensity of the impact relationship between on the ESs, namely, land urbanization > population urbanization > economic urbanization. In 2018, this trend was reversed, with economic urbanization (GDPD, Moran’s I = −0.407) and population urbanization (POP, Moran’s I = −0.549) experiencing a significant increase and a small decrease, respectively, in their influence on the ESs, altering the ranking of the strength of the influence relationship on the ESs (land urbanization > population urbanization > economic urbanization).

3.3.2. Local Bivariate Moran’s I Analysis

To explore the spatial response of the ESs to urbanization disturbances in the study area, we also carried out local correlation analysis of the variables. The local bivariate LISA results revealed four spatial correlations between urbanization and ESs. The results of spatial agglomeration mapping indicated that the effects of population growth, GDP growth and construction land expansion on the ESs in the study area were similar and revealed a slight change in 2000/2005/2015 and a large change in 2018.

The spatial correlation between the construction land share and ESs is shown in Figure 7, with L-L type agglomeration areas mainly occurring in Beichen District, Xiqing District and Jinan District of Tianjin in 2000 and the Bashang Plateau Region from 2000–2015 and disappearing in 2018. The H-L agglomeration areas in the study area were concentrated in the central urban areas of Beijing and Tianjin in the study years, where the H-L agglomeration areas in Beijing did not change and those in Tianjin expanded from the central urban area toward the periphery with increasing urban expansion. The L-H concentration areas in the study area were concentrated in the region of the Yanshan Mountains.

The spatial correlation between the GDP density and ESs is shown in the figure below, with L-L type agglomeration areas largely concentrated in the Bashang Plateau Region from 2000 to 2015, L-L type agglomeration areas emerging near the central city of Tianjin in 2015, L-L type agglomeration areas disappearing in the Bashang Plateau Region in the study area in 2018, and L-L type agglomeration areas expanding outwards along the periphery of the central city of Tianjin. The H-L type aggregation areas in the study area were concentrated in the central urban areas of Beijing and Tianjin in the study years, with the H-L type aggregation area in Beijing changing non-significantly and that along the periphery of the central urban area of Tianjin beginning to change from an H-L type aggregation area to an L-L type aggregation area.

The spatial correlation between the population density and ESs is shown in the figure below. During the study period, the H-L type agglomeration areas in the study area were concentrated in the central cities of Beijing and Tianjin, where the Beijing H-L type agglomeration area did not change, while from 2000–2015, those in Tianjin city and the four districts of Dongli, Xiqing, Jinan and Beichen experienced transformation into H-L and L-L type agglomeration areas, respectively. Finally, in 2018, this area and the Binhai New Area all became L-L type agglomeration area. In addition, L-L type aggregation areas were also present in large numbers in the Bashang Plateau Region from 2000–2015 but disappeared in this region in 2018.

3.3.3. Spatially Nonstationary Effects of Urbanization on the ESs

To further explore the spatial response of the ESs to each type of urbanization indicator in each district and county, this paper analyzed the spatial heterogeneity in the influence of the different urbanization factors on the ES in each district and county within the BTH separately with the GWR model and visualized the regression coefficients of each influencing factor with the quantile classification method to further assess the differences in the spatial responses of the ESs to urbanization in the different districts and counties. The final parameters of the GWR model are listed in

Table 3. Calculation results of the geographical weighted regression (GWR) model.

Years	2000			2005			2015			2018
Variables	CLP	GDPD	POPD	CLP	GDPD	POPD	CLP	GDPD	POPD	CLP	GDPD	POPD
R square	0.79	0.83	0.71	0.85	0.84	0.81	0.83	0.55	0.71	0.84	0.63	0.71
Adjusted R square	0.76	0.81	0.67	0.83	0.82	0.78	0.81	0.49	0.67	0.81	0.58	0.67
DIFF of Criterion	−17.10	−60.22	−27.54	−15.64	−43.72	−29.48	−12.86	−17.11	−14.61	−17.70	−4.76	−14.40

. R² indicates the model fit: high values indicate a high model performance. The DIFF of Criterion was used to test for the presence of spatial non-stationarity. If the DIFF of Criterion value is negative, there then exists spatial variability.

The results of the ESs responses to each urbanization indicator during the study period and the amount of change in the intensity of each urbanization indicator driving the ESs are shown in Figure 8 and Figure 9, respectively. The proportion of the built-up area in each district and county was significantly and negatively correlated with the ESs in the study area. A high intensity of the driving effect of the construction land area on regional ESs was concentrated in areas with strong ESs in the Yan Mountains and relatively low levels of spatial urbanization. In contrast, the regression coefficients were lower in the plain regions with better ESs and a higher proportion of urban construction areas and in the Bashang Plateau Region with low coefficient of the ESs and a smaller proportion of built-up areas, where the strength of the driving effect of spatial urbanization on the ESs was lower. The negative driving intensity of CLP on ESs showed four different characteristics of “enhanced-weakened-enhanced”, “weakened-enhanced-weakened”, “weakened-enhanced” and “enhanced-weakened” in Yanshan-Taihang Mountain region, Bashang Plateau region and plain areas, respectively. The incremental GDP density was negatively correlated with regional ESs. During the study period, the region with a stronger ESs driving force the GDP density was concentrated in Beijing and its northern districts and counties in 2000, which shifted to the Tangshan region in 2005 and gradually decreased to the west, and this region then shifted to the southern region of Baoding city in 2015. However, the spatial pattern of the regional driving intensity did not change significantly in 2018, but the driving intensity weakened in the Yanshan-Taihang Mountains region and the Bashang Plateau region, and strengthened in the plains region. The population density was negatively correlated with the ESs in the study area. During the study period, the intensity of the impact of population density on ESs is characterized by an “increasing-weakening-increasing” pattern in the study area. Driving intensity changed from 2000 to 2005, with the center of the high-driving intensity region shifting from the Yan and Taihang Mountains to Tianjin and surrounding local areas and decreasing outwards. From 2005–2018, the regional center of high coefficient shifted back to the Yan and Taihang Mountains region, exhibiting spatial distribution characteristics of high coefficient in the mountains and low coefficient in the plains.

4. Discussion

4.1. ES Spatiotemporal Evolution

In this study, the critical ESs were quantified using the integrated ES index in 2000, 2005, 2015 and 2018., which is more accurate than the results of the best-effort land use type (ecosystem service type) assessment system using the equivalence factor method, reflecting the richness of information inherent in ESs at the spatial and temporal scales due to natural factors [68,69]. During the course of overall ESs change in the BTH, ESs first improved during the 2000–2015 period and declined during the 2015–2018 period. As the political center of China, regional development is closely related to national policies. After winning the right to host the 2008 Summer Olympics in Beijing in 2001, the Chinese government implemented measures to shut down heavily polluting factories surrounding Beijing, while large-scale ecological restoration projects such as returning farmland to forestland and grassland, as well as the construction of regional greening and recreational support facilities, improved Ess [70], resulting in a significant improvement from 2000–2015. In 2015, the Chinese government adopted the Beijing-Tianjin-Hebei Synergistic Development Planning Outline, which promoted the synergistic development of the BTH as a major national strategy, which undoubtedly accelerated the development of the BTH. Although the government actively carried out construction of the Beijing-Tianjin-Hebei Ecological Coordination Circle during this period, this was accompanied by rapid urbanization, which led to a gradual decline in the intensity of regional ESs. In the spatial dimension, high-ES areas were concentrated in the Yan Mountain Range region due to the mountainous environment limiting human activities [71], while low-ES areas were concentrated in the Bashang Plateau Region because this region is located at the center of the northern agropastoral interlacing zone and is a fragile ecological zone in northern China [72,73], which faces environmental problems such as land degradation and soil erosion with increasing human activities. The plain areas fluctuated more notably during the study period, and areas adjacent to urban centers attained a lower intensity of their ESs and were more vulnerable to anthropogenic impacts [74].

4.2. Spatial Heterogeneity Analysis of the Impact of Urbanization on ESs in BTH Urban Agglomerations

Human activities are considered a major driver of global environmental change, and rapid regional urbanization in particular exerts a significant impact on the health of ecosystems themselves and their sustainable development [75,76]. In this study, the impact of overall regional urbanization on ESs was not only explored, but the impact of urbanization in each district and county within the selected urban agglomeration on its own and surrounding region-wide ESs intensity was also explored and analyzed in depth. The results of global bivariate Moran’s I indicate that the GDP density, population density and proportion of land used for construction all exerted a negative impact on ESs, while the negative impact on the county ESs was greater than that of the GDP density, population density, and proportion of land used for construction. This result could occur because the urbanization process is accompanied by population concentration, economic development and expansion of built-up land and because urban built-up land expansion often more directly impacts ESs by directly altering the nature of land. To gain a more comprehensive understanding of the effects of urbanization on the ESs, this paper also measured the spatial correlation between urbanization and ESs in neighboring areas through global bivariate Moran’s I, and the results reveal that the high-urbanization value areas in the central cities of Beijing and Tianjin exerted a negative impact on the ESs in neighboring areas, with a significant spillover effect. The reason for this phenomenon may be that Beijing and Tianjin require large amounts of energy and materials for the twin engines of development of the BTH, and the neighboring areas, as the energy and material suppliers of these highly urbanized areas, continuously provide raw materials for urban development by exploiting natural resources and emitting pollutants including harmful gases such as SO_2, garbage and industrial waste, which generates a negative impact on the surrounding ecosystem [77,78,79]. This is why we observed large areas of low-ES surrounding large, highly concentrated urban centers.

Compared with Geographical detector application and ordinary least squares, the regression using GWR model is more beneficial to consider and can reflect the intensity of urbanization driving ESs more intuitively [80]. The GWR results demonstrate that the stronger driving force of CLP on the ESs during the study period was concentrated in the Yan and Taihang Mountain Range region with high-ES coefficient, mainly because at the early stage of construction, the urban built-up area tends to expand radially, and rapid expansion of construction land brings about the destruction of ecosystems [81], which could be more significant in areas with better ecosystems. The driving force of CLP in cities at the middle and late stages of urban development is much weaker. In 2005 and 2015, there occurred clear spatial changes in the intensity of ESs enhancement per unit area of the GDP density and POP density. The change in 2005 was attributed to the beginning of the industrial evacuation process of Beijing toward neighboring cities, of which Tianjin and Tangshan, as an important part of the Beijing-Tianjin-Tangshan industrial zone and with a more solid industrial base, received a large number of polluting factories such as iron and steel plants and experienced a notable population influx. This could undoubtedly cause disturbance and damage to the regional ecosystem [82,83]. New development of Baoding as part of the BTH has attracted a large number of enterprises and employees, causing a very high pressure on the ESs. Moreover, as the industrial model of the BTH is gradually transformed from secondary to tertiary industries, GDP enhancement will rely more on the development of new industries, high-tech development and services industries, with a lower spatial demand and environmental pressure, and the driving ability of the GDP density on the ESs will decrease [84,85].

4.3. Strategies for Optimizing Ecological Resources and Spatial Control

To mitigate the negative impacts of BTH urbanization on the ESs, we should develop strategies to maintain and improve the capacity of the ecosystem so that the ESs along a positive direction. Based on the results of this study, we propose the following control strategies: (1) For the two megaregions of Beijing and Tianjin, their high value areas of urbanization in the study area and the critical ESs that maintain regional ecological security are low, which will undoubtedly make the region’s face higher security risks. In the future, the impact on the surrounding ecosystems should be reduced by decongesting nonessential urban functions, controlling population clustering and urban expansion in this region. Meanwhile, green areas and greenways should be strengthened in areas with high CLP values to enhance the capacity of ESs. (2) In areas with high-ES, such as the Yan Mountains, it is important to enhance the protection of the ecosystems in the area. Meanwhile, land and population urbanization in this area have a strong driving effect on ESs, so the urbanization of land and population should be strictly controlled and managed to prevent the rapid depletion of natural resources in this area [86]. This is not only a matter concerning its own ecological environment but also a matter concerning the protective barrier for ecological security and healthy development of Beijing, the capital of China. Therefore, we should improve the cooperation mechanism for joint prevention and control of regional pollution, actively implement ecological compensation measures, and enhance technical and financial support for ecological protection in the region [87,88]. (3) In the Bashang Plateau region and some plains urban cities, which have a low-ES and are also less disturbed by urbanization, the sustainability of the ecosystem should be emphasized in the future, and the regional ESs should be incorporated into the norms of sustainable urban development and management, and strictly control the development of polluting industries and large-scale urbanization construction to prevent them from causing damage to regional ecosystems [89,90].

4.4. Research Deficiencies and Prospects

Despite the essential information provided, this study suffers certain limitations on a limited data precision in this paper. The data precision we obtained in this study is 1 km × 1 km. A coarse precision may lead to errors in the calculated spatial data of the ESs, and limited by the accessibility to urbanization-related data, regression analysis was conducted only at the district and county level. Therefore, we suggest obtaining data with a higher precision to analyze the spatial response of ESs to different types of urbanization indicators on a fine spatial scale in the future. Meanwhile, the relationship between ESs and urbanization is complex, and we hope to add more qualitative indicators reflecting human well-being and better and more comprehensive model analysis the relationship between urbanization and ESs in depth.

5. Conclusions

Because of accelerated urbanization, human activities have greatly changed the ESs of urban agglomerations. Analyzing the evolutionary characteristics of ESs and the impact of urbanization on the ecosystem supply and demand is essential for regional ecological management and sustainable development. This study was based on the BTH in China. A long time span of nearly 20 years from 2000–2018 was chosen to analyze the spatiotemporal evolution of ecological services in the BTH, and on this basis, the bivariate global Moran’s I and bivariate local Moran’s I were applied to explore the spatial correlation between different types of urbanization indicators and ESs and the location where the spatial aggregation phenomenon occurred. Finally, the GWR model was employed to examine the spatial ESs response to different types of urbanization indicators.

The results were as follows: first, the ecological services in the BTH exhibited an increasing trend and a subsequent decreasing trend from 2000 to 2018, reaching its highest index in 2015. Second, the ESs from 2000 to 2018 indicated obvious spatial divergence characteristics with the Yan and Taihang Mountains region as boundaries, a weakening trend in the entire northwestern region, and a gradual decrease in the southeastern region since 2015. In the process of exploring the spatial response of ESs to urbanization with the bivariate global and local Moran’s I, it was found that the proportion of built-up land, POP density and GDP density all exhibited relatively obvious negative driving characteristics, with the proportion of built-up land generating the greatest negative driving force on the ESs. Moreover, in terms of the spatial distribution, the areas around the central cities of Beijing and Tianjin, the Bashang Plateau Region and the Yan Mountains already revealed a clear aggregation effect. Fourth, the GWR model results suggest that the proportion of built-up land imposed a greater influence on areas with high-ES, while the spatial response of the population density and the intensity of the effect of the GDP density were largely influenced by notable and highly polluting, energy-intensive industries. Based on the analysis results, the mechanisms of the urbanization impact on ESs were explored, and future policy recommendations were formulated to mitigate the impacts of urban expansion on ESs. This is of value to the BTH in achieving suitable ecosystem management.