Vegetation Dynamics and Their Response to the Urbanization of the Beijing–Tianjin–Hebei Region, China

: Rapid global urbanization has caused substantial changes in land cover and vegetation growth. Rapid urban growth in a short time has escalated the conﬂicts between economic development and ecological conservation, particularly in some metropolitan regions. However, the e ﬀ ects of rapid urbanization on vegetation have not been fully captured, especially accounting for the latest ecological development initiatives. In this study, we chose a typical urban agglomeration, the Beijing–Tianjin–Hebei (BTH) urban agglomeration in China, and analyzed the vegetation variation and the impacts of urbanization on the vegetation growth based on transferable methods, using data such as the Normalized Di ﬀ erence Vegetation Index (NDVI) and the nighttime light (NTL). The results indicate signiﬁcantly enhanced vegetation growth in the BTH region, with a strikingly spatial pattern of greening in the northwest, and browning in the southeast from 2001 to 2018. Besides this, the results enclose most of the areas (72%) of built-up land in the BTH, which tended to brown in the process of rapid urban development, while 27% greened with increasing urbanization. This means that the vegetation’s response to urbanization shows apparent di ﬀ erences and geographic heterogeneity along the urbanization gradient at the urban agglomeration scale. Parts of the periphery of the metropolis and the central areas of developing cities may experience a browning trend; however, the core urban areas of urbanized metropolises demonstrate greening, rather than browning. Furthermore, this study provides solid evidence on the remarkable greening impacts of several ecological restoration projects which are currently underway, especially in ecologically fragile areas (e.g., the suburbs). The implications derived from the urban ecological development and the transferable methodology deployed in this paper facilitate the unfolding relationships between urbanization and social-ecological development. Our ﬁndings provide new insights into the interactions between vegetation dynamics and urbanization at the regional level.


Introduction
Urbanization is a common phenomenon in the process of global modernization, and it is an important indicator of social progress and economic development [1]. Urban development usually leads to dramatic land use alternation, and it is usually recognized as a process that can cause ecological pressures to terrestrial ecosystems [2,3] resulting in irreversible consequences [4][5][6], such as environmental pollution, the shortage of resources, a drastic reduction of wetlands, the loss of dynamics and explore the inner mechanism of the effects of human activities on vegetation changes at the agglomeration scale on the basis of continuing urban morphological observations and vegetation evolution trends.
As China's third economic growth pole right after the Yangtze River Delta urban agglomeration and the Pearl River Delta urban agglomeration, the Beijing-Tianjin-Hebei (BTH) urban agglomeration-the 'Capital Economy Circle' of China-is a typical region in which to undertake this investigation in order to ascertain the extent to which human activities have influenced urban vegetation. Moreover, the BTH urban agglomeration has developed dramatically into a national political, economic, and cultural center since the 1990s [42,43]. However, the continuous expansion of its urban agglomerations, especially Beijing and Tianjin, has led to complex problems in local ecosystems [44], which not only pose a threat to ecological security but may also substantially hinder regional sustainability development [45]. In fact, in the progress of development and the transformation of industry, it is widely accepted that great attention should be paid to the promotion of an excellent ecological environment, developing a green economy, and building a 'Two-Oriented Society' (i.e., a resource-saving and environment-friendly society). The restoration of ecological systems and the economic development of the BTH region should be addressed simultaneously in the long term. Nevertheless, at present, few studies have addressed this issue at the same time from the point of view of the regional synthetic sustainable development. This paper aims to fill this gap by using remote sensing techniques. Accordingly, based on the analysis mentioned in the above literature, this study proposes a feasible approach to measure vegetation and urbanization variation in the BTH. This study has three components. First, using the Moderate Resolution Imaging Spectroradiometer (MODIS) image, we analyzed the spatiotemporal evolution of the vegetation coverage in BTH from 2001 to 2018. Second, based on the NTL data, the urbanization evolution trend of the region from 2001 to 2013 was analyzed. Finally, we utilized the spatial division of the vegetation variation trend combined with nighttime light intensity variation trend to identify the internal mechanism of the vegetation dynamics and urbanization evolution in the different regions.
The primary purpose of the study is to explore the characteristics of the vegetation's response in urban agglomerations, and to investigate the relationship between urban development and vegetation dynamics in cities with different development intensities belonging to the same urban agglomerations. This study can offer a scientific basis for the sustainable development of urban agglomerations with similar features all over the world.

Study Area
The BTH region-the Capital Economy Circle of China-lies between 36 • 03 -42 • 40 N and 113 • 27 -119 • 50 E in the North China Plain, with a total area of 218,000 km 2 . It includes two municipalities directly under the central government (i.e., Beijing and Tianjin) and 11 cities at prefecture level in Hebei Province ( Figure 1). The elevation is higher to the northwest and lower to the southeast (Figure 1), and there are three major geomorphic units: Zhangjiakou Bashang Plateau, the Yanshan-Taihang Mountains, and Hebei Plain. The typical temperate monsoon climate in this region is characterized by high temperatures and rain in the summer, and coldness and dryness in the winter. There are various types of vegetation in the BTH region. Herbs dominate the grassland area in the northern plateau, and the main vegetation in the Yanshan-Taihang Mountains is shrubs and bushes. The vegetation types in the mountain basins in northwestern Hebei are mainly forests and shrubs. Most of the southern plains are cultivated land and artificial urban vegetation.

Vegetation Index Data
The remote sensing data used in this study include the MODIS NDVI dataset and the DMSP/OLS dataset. We used the Application for Extracting and Exploring Analysis Ready Samples (AppEEARS) (https://lpdaacsvc.cr.usgs.gov/appeears/) to download the BTH urban agglomeration MOD13Q1 NDVI products from 2001 to 2018, which is a 16 day synthesis with 23 images per year, within a spatial resolution of 250m. This product has not only eliminated the impact of sensor performance attenuation [46] but also improved the integration algorithm of the vegetation index products, and the product data quality is higher than the previous versions. The data downloaded through AppEEARS were consolidated into the WGS84 coordinate system. In order to eliminate the noise phenomenon caused by clouds, snow, and adverse atmospheric conditions, the maximum synthesis method (MVC) was further used to generate monthly datasets. In order to reflect the vegetation information more accurately, we limited the NDVI observation period to the Chinese major plant growth season (i.e., April to October), and used it to calculate the annual NDVI value.

Nighttime Light Data
We acquired the DMSP/OLS (Version 4) stable light annual image dataset (2001-2013) from the National Oceanic and Atmospheric Administrationʹs National Geophysical Data Center (NOAA/NGDC) website (http://www.ngdc.noaa.gov/). The stable light image is an annual raster image containing information about the average light intensity at night. The images include steady light sources in cities, towns, and other places, and remove the influence of unstable noise, such as moonlight clouds, light, and fire. These data products were captured from four different satellites: F14, F15, F16, and F18. They were generally acquired at 9-10 p.m. for each local time, with a spatial resolution of 30 radians (about 1 km).
Since the resolution of the NTL data is 30 radians, the image grid decreases as the latitude increases [47]. Thus, we converted the projection coordinate systems of all of the images into the Lambert projection coordinate system, in order to avoid the deviation of the image grid deformation. There are differences in the orbit parameters, sensing equipment, and atmospheric refraction between the various satellites, so the NTL values acquired by the different satellites are not consistent [48]. Consequently, it is essential to preprocess the data in order to mitigate the negative impacts, such as

Vegetation Index Data
The remote sensing data used in this study include the MODIS NDVI dataset and the DMSP/OLS dataset. We used the Application for Extracting and Exploring Analysis Ready Samples (AppEEARS) (https://lpdaacsvc.cr.usgs.gov/appeears/) to download the BTH urban agglomeration MOD13Q1 NDVI products from 2001 to 2018, which is a 16 day synthesis with 23 images per year, within a spatial resolution of 250 m. This product has not only eliminated the impact of sensor performance attenuation [46] but also improved the integration algorithm of the vegetation index products, and the product data quality is higher than the previous versions. The data downloaded through AppEEARS were consolidated into the WGS84 coordinate system. In order to eliminate the noise phenomenon caused by clouds, snow, and adverse atmospheric conditions, the maximum synthesis method (MVC) was further used to generate monthly datasets. In order to reflect the vegetation information more accurately, we limited the NDVI observation period to the Chinese major plant growth season (i.e., April to October), and used it to calculate the annual NDVI value.

Nighttime Light Data
We acquired the DMSP/OLS (Version 4) stable light annual image dataset (2001-2013) from the National Oceanic and Atmospheric Administration's National Geophysical Data Center (NOAA/NGDC) website (http://www.ngdc.noaa.gov/). The stable light image is an annual raster image containing information about the average light intensity at night. The images include steady light sources in cities, towns, and other places, and remove the influence of unstable noise, such as moonlight clouds, light, and fire. These data products were captured from four different satellites: F14, F15, F16, and F18. They were generally acquired at 9-10 p.m. for each local time, with a spatial resolution of 30 radians (about 1 km).
Since the resolution of the NTL data is 30 radians, the image grid decreases as the latitude increases [47]. Thus, we converted the projection coordinate systems of all of the images into the Lambert projection coordinate system, in order to avoid the deviation of the image grid deformation. There are differences in the orbit parameters, sensing equipment, and atmospheric refraction between Sustainability 2020, 12, 8550 5 of 21 the various satellites, so the NTL values acquired by the different satellites are not consistent [48]. Consequently, it is essential to preprocess the data in order to mitigate the negative impacts, such as unavoidable noises, created during the data's acquisition. In this study, we integrated remote sensing data from multiple satellites and used the NTL data systematic correction method proposed by Cao et al. to correct the data [49], in order to improve the consistency and comparability of the NTL data for our specific research time.

Additional Data For the Research Area
The vector data of ecological engineering boundary in the BTH region, the land cover data, and the vector data of the administrative boundary of the BTH urban agglomeration were obtained from the Resource and Environment Data Cloud Platform of the Chinese Academy of Sciences (http://www.resdc.cn/DataSearch.aspx). The statistical data and the green coverage of the urban built district were collected from the China City Construction Statistical Yearbook (2002-2018) (http://www.mohurd.gov.cn/).

The Calculation of the Fractional Vegetation Coverage
There is a linear relationship between the NDVI and the fractional vegetation coverage (FVC) [50]. Thus, the vegetation coverage information can be calculated by establishing its conversion relationship. In this study, we used a remote sensing estimation method based on a pixel dichotomous model design, which set NDVI as a parameter in order to calculate the vegetation coverage. It does not need to estimate the leaf area index (LAI) and other parameters that require complicated derivation, and is suitable for different vegetation types. The FVC was defined as follows: where NDVI soil indicates the bare land, and is equivalent to the minimum value; and NDVI veg is the NDVI value of a completely vegetation-covered pixel, which represents the maximum value. In fact, due to the influence of weather conditions, vegetation types, seasonal changes, and other factors, NDVI soil and NDVI veg may fluctuate in different years. The preprocessed NDVI data was used as the mixed pixel vegetation index value. We set a 5% confidence level as the upper and lower thresholds of the monthly average NDVI index to capture the NDVI veg and NDVI soil value as previous studies [50].

FVC and NTL Trend Detection
We took advantage of one-dimensional linear regression analysis to study the variation trend of the FVC quantitatively. This method can illustrate the variation trend of each grid, and thus comprehensively reflect the regional spatial and temporal heterogeneity. The trend slope of each pixel was calculated by the ordinary least square (OLS) method: (2) in which n means the cumulative number of years, FVC i is equal to the FVC value of the i-th year, and slope indicates the variation trend of the FVC values during the observation period. Slope > 0 means that the vegetation tends to develop in a positive trend during the monitoring period; otherwise, slope < 0 means that the vegetation growth is destroyed, and that the vegetation status tends to Sustainability 2020, 12, 8550 6 of 21 deteriorate [51]. Additionally, we tested the significance of the trends for all of the pixels by using the F test as follows: (y i −ŷ i ) 2 is the regression squared sum, y i is the pixel value of the i-th year,ŷ i is the regression value, y is the average value of coverage during the monitoring period, and n is the number of monitoring years. According to the test results, we divided the changes into five levels: Significant greening (slope > 0, p ≤ 0.05); Green (slope > 0, p > 0.05); Significant browning (slope < 0, p ≤ 0.05); Browning (slope < 0, p > 0.05); and No change (slope = 0). In addition, we utilized the former formulas to calculate the variation trend of the NTL in the BTH urban agglomeration. Since most areas showed a growing trend, we used the Jenks Natural Breaks Function to group the positive values to distinguish the degree of urbanization in the different zones (i.e., dramatic increase, moderate increase, slight increase, no change, and decrease).

Calculation of the Average Nighttime Light Intensity Index
The level of urbanization is a comprehensive reflection of many factors (e.g., population, industrial structure, geographical space). Similarly, the DMSP/OLS nighttime light index is a comprehensive reflection of the interaction of these factors. In this study, the average nighttime light intensity index (ANLI) was used as a characteristic indicator of the urbanization level, in order to identify the urbanization level of 13 cities in the BTH region. The ANLI was calculated as follows: (4) in which ANLI represents the average nighttime light intensity index; DN i is the light intensity of the i-th level in the area; n i is the total number of pixels of the i-level light intensity in the area; N indicates the number of all of the light pixels (1 ≤ N ≤ 63) in the area; and N max is the maximum DN value.

Correlation Analysis of the Urbanization Level and Vegetation Dynamics
We established a correlation model between the DN value and FVC in order to analyze the relationship between the urbanization level and the vegetation coverage at the administrative district and pixel levels, respectively. At the pixel level, the same interval division method was applied. This approach has been proven to mitigate the influence of noise pixels on the mathematical results [52]. We subdivided the DN value into six layers at intervals of 10 (i.e., 1-10, 20-30, . . . , 50-63), and then randomly selected 50 sampling points from each sample interval in order to promote the uniformity of sampling, and recorded the corresponding DN value and FVC value of each sample pixel. After the sampling, we implemented a regression model and a significance test to evaluate whether the relevance was significant between the urbanization level and the vegetation coverage. If the correlation coefficient is negative, it means that the DN value has a negative correlation with the vegetation coverage of the BTH region; otherwise, DN has a positive effect on FVC. Additionally, in order to show the spatial heterogeneity of the coupling result, both the ArcGIS and Matlab platforms were applied in order to execute the spatial superposition and statistical analysis of the NTL and FVC in the BTH region from 2001 to 2013. This method can illustrate the relationship between urbanization and vegetation variations in various contexts.

Vegetation Variation Analysis of the BTH Region
Based on the NDVI time-series data from 2001 to 2018, the spatial characteristic of the average vegetation coverage of the BTH urban agglomeration was calculated ( Figure 2a). Meanwhile, the average FVC of each city from 2001 to 2018 was counted, and the significant differences amongst the cities were observed, as shown in Figure 2b.

Vegetation Variation Analysis of the BTH Region
Based on the NDVI time-series data from 2001 to 2018, the spatial characteristic of the average vegetation coverage of the BTH urban agglomeration was calculated ( Figure 2a). Meanwhile, the average FVC of each city from 2001 to 2018 was counted, and the significant differences amongst the cities were observed, as shown in Figure 2b. As shown in Figure 2a (the proportion of the vegetation coverage), the vegetation coverage in the BTH region was fairly good. Areas with an average FVC of less than 40% only account for 8.18% of the total area, while the areas with 0.5 < FVC < 0.7 account for 53.18%, and the areas with high vegetation coverage (FVC > 0.7) account for 25.41%. Overall, areas with vegetation coverage greater than 50% account for 78.59%. According to Figure 2a,b, the vegetation coverage within the BTH region presents significant spatial heterogeneity, especially in the western and eastern parts of the BTH region. Therefore, the results, shown in Figure 2, indicate that there was severe land degradation in Zhangjiakou, which is located in the ecologically fragile area. Moreover, there were abundant mudflats in Tianjin and Cangzhou, resulting in their FVC being below the average level. The FVC of the plains in the southeast of the BTH region was moderate. In this region, the low vegetation coverage areas were mainly distributed among the built-up areas, which had suffered from intensive human activities.  As shown in Figure 2a (the proportion of the vegetation coverage), the vegetation coverage in the BTH region was fairly good. Areas with an average FVC of less than 40% only account for 8.18% of the total area, while the areas with 0.5 < FVC < 0.7 account for 53.18%, and the areas with high vegetation coverage (FVC > 0.7) account for 25.41%. Overall, areas with vegetation coverage greater than 50% account for 78.59%. According to Figure 2a,b, the vegetation coverage within the BTH region presents significant spatial heterogeneity, especially in the western and eastern parts of the BTH region. Therefore, the results, shown in Figure 2, indicate that there was severe land degradation in Zhangjiakou, which is located in the ecologically fragile area. Moreover, there were abundant mudflats in Tianjin and Cangzhou, resulting in their FVC being below the average level. The FVC of the plains in the southeast of the BTH region was moderate. In this region, the low vegetation coverage areas were mainly distributed among the built-up areas, which had suffered from intensive human activities. According to the slope value of the FVC in the BTH region from 2001 to 2018 (Figure 3a), the ratio of various pixels was counted (Figure 3c). In the past 18 years, the vegetation coverage in most areas of BTH varied to some extent; the area with no change only accounts for 0.99%; the area with an increasing FVC accounted for 55.13%, in which the area with significant restoration (p ≤ 0.05) accounts for 18

Analysis of the Nighttime Light Characteristics
We modified the DMSP/OLS nighttime light data from 2001 to 2013 according to the previous correction method, and obtained 13 NTL maps of BTH ( Figure 4). On the whole, the nighttime light intensity value was always high in the southeast and low in the northwest from 2001 to 2013. In this period, the nighttime light index apparently expanded, especially in Beijing and Tianjin. Simultaneously, the NTL can be utilized for morphological observations of built-up districts in the region. The built-up areas in the southeast were densely clustered, while those in the northwest were sporadic.

Analysis of the Nighttime Light Characteristics
We modified the DMSP/OLS nighttime light data from 2001 to 2013 according to the previous correction method, and obtained 13 NTL maps of BTH (Figure 4). On the whole, the nighttime light intensity value was always high in the southeast and low in the northwest from 2001 to 2013. In this period, the nighttime light index apparently expanded, especially in Beijing and Tianjin. Simultaneously, the NTL can be utilized for morphological observations of built-up districts in the region. The built-up areas in the southeast were densely clustered, while those in the northwest were sporadic.
In order to compare the development degree of different cities clearly, and to better observe the outstanding clusters in Figure 3, we obtained the average night light density index (ANLI) heat map of each city for 13 years ( Figure 5). The figure demonstrates that the average light density index of each city increased steadily in the 13 years, but there was an apparent difference among 13 cities. Beijing and Tianjin showed distinct development of their ANLI in the year of 2013, which reached 0.56 and 0.47, respectively. Besides this, due to the radiation effects of the two metropolises (i.e., Beijing and Tianjin), Langfang had the highest level of ANLI compared to the other cities in the Hebei province.
Although the values of the ANLI index of Hengshui, Chengde, and Cangzhou rose continuously, they were still at the lowest level within the BTH region. In order to compare the development degree of different cities clearly, and to better observe the outstanding clusters in Figure 3, we obtained the average night light density index (ANLI) heat map of each city for 13 years ( Figure 5). The figure demonstrates that the average light density index of each city increased steadily in the 13 years, but there was an apparent difference among 13 cities. Beijing and Tianjin showed distinct development of their ANLI in the year of 2013, which reached 0.56 and 0.47, respectively. Besides this, due to the radiation effects of the two metropolises (i.e., Beijing and Tianjin), Langfang had the highest level of ANLI compared to the other cities in the Hebei province. Although the values of the ANLI index of Hengshui, Chengde, and Cangzhou rose continuously, they were still at the lowest level within the BTH region. By analyzing the slope trend and verifying its significance (p < 0.05), we obtained the distribution map of the NTL variation slope in the BTH urban agglomeration. Based on the above research results, we recognized Hengshui City-which had the lowest ANLI index in the BTH region-as the sample city, and selected the built-up area statistics data closest to the NTL data as the optimal threshold to extract the city built-up areas range ( Figure 6). Overall, the NTL increased significantly in the southeast of BTH. However, in the northern and northwest regions, only Chengde and Zhangjiakou had a particular increase in their urban areas. During the past 13 years, the areas with the most By analyzing the slope trend and verifying its significance (p < 0.05), we obtained the distribution map of the NTL variation slope in the BTH urban agglomeration. Based on the above research results, we recognized Hengshui City-which had the lowest ANLI index in the BTH region-as the sample city, and selected the built-up area statistics data closest to the NTL data as the optimal threshold to extract the city built-up areas range ( Figure 6). Overall, the NTL increased significantly in the southeast of BTH. However, in the northern and northwest regions, only Chengde and Zhangjiakou had a particular increase in their urban areas. During the past 13 years, the areas with the most energetic light intensity growth were mainly the peripheral areas of Beijing and Tianjin. Meanwhile, Tangshan, Shijiazhuang, Handan, and Xingtai also experienced a significant growth of their peripheral light intensity. The cluster expansion showed that the connectivity between the main urban areas and the surrounding counties dramatically increased. Nevertheless, the other cities mainly showed a dotted growth pattern of NTL in their urban areas and counties, and the connectivity between the urban areas and counties was not apparent.

Correlation Analysis of NTL and FVC
The regression analysis results of the mean DN and FVC are shown in Figure 7. The regression results of the mean DN and FVC of the 13 cities and 300 sample points both demonstrate that the FVC had a negative correlation with the DN value in the BTH region, and the correlation coefficients were respectively 0.404 and 0.675. The results mean that areas with high NTL generally have low vegetation coverage. Because the NTL is a reflection of the urbanization degree in human social and economic activities, regional urbanization has a significant spatial correlation with vegetation degradation. That is to say, a higher intensity of human activities may worsen the condition of the vegetation.

Correlation Analysis of NTL and FVC
The regression analysis results of the mean DN and FVC are shown in Figure 7. The regression results of the mean DN and FVC of the 13 cities and 300 sample points both demonstrate that the FVC had a negative correlation with the DN value in the BTH region, and the correlation coefficients were respectively 0.404 and 0.675. The results mean that areas with high NTL generally have low vegetation coverage. Because the NTL is a reflection of the urbanization degree in human social and economic activities, regional urbanization has a significant spatial correlation with vegetation degradation. That is to say, a higher intensity of human activities may worsen the condition of the vegetation. In order to make the time interval consistent, we calculated the FVC variation slope from 2001 to 2013. We overlaid the FVC slope and the NTL slope results of BTH to reveal the coupling relationship between the trends of the FVC and NTL of the BTH urban agglomeration. The results are shown in Figure 8, and are divided into four categories: Category I indicates that the vegetation coverage increased in the process of urban development; Category II indicates that the vegetation coverage increased to a certain extent while the city shrank; Category III refers to the destruction of vegetation and the decline of the vegetation coverage in the process of counter-urbanization; Category IV refers to the decrease of the vegetation coverage in the process of urbanization, considering that the area initially used for plant growth was occupied by construction land. According to the above classification, the grid data were statistically analyzed at the pixel scale ( Figure 9). Due to the influence of outliers, the original statistical results made it difficult to show the partition. We thus excluded the exceptional values of two variables in order to show the distribution of the various pixels more intuitively.  In order to make the time interval consistent, we calculated the FVC variation slope from 2001 to 2013. We overlaid the FVC slope and the NTL slope results of BTH to reveal the coupling relationship between the trends of the FVC and NTL of the BTH urban agglomeration. The results are shown in Figure 8, and are divided into four categories: Category I indicates that the vegetation coverage increased in the process of urban development; Category II indicates that the vegetation coverage increased to a certain extent while the city shrank; Category III refers to the destruction of vegetation and the decline of the vegetation coverage in the process of counter-urbanization; Category IV refers to the decrease of the vegetation coverage in the process of urbanization, considering that the area initially used for plant growth was occupied by construction land. According to the above classification, the grid data were statistically analyzed at the pixel scale ( Figure 9). Due to the influence of outliers, the original statistical results made it difficult to show the partition. We thus excluded the exceptional values of two variables in order to show the distribution of the various pixels more intuitively.
Among the four categories, the proportion of Category IV was the highest, with a total area of 23,054.42 km 2 , accounting for 72.27% of the built-up area of the BTH region; the proportion of Category I was the second, accounting for 26.89%, and 8281.73 km 2 ; the proportions of Category II (0.34%) and Category III (0.50%) were extremely tiny. This shows that the urbanization level in most of the built-up areas of BTH increased in the past 18 years, but there was apparent spatial heterogeneity in the vegetation growth in the different regions. With regard to Category I, it has two distinct distribution forms. On the one hand, it is distributed in the core areas of metropolises (e.g., Beijing, Tianjin, and Shijiazhuang). On the other hand, Category I was also located in the suburbs or coastal areas of Zhangjiakou, Chengde, Tianjin, and Hengshui. Category IV was mainly distributed in suburban areas, such as Beijing, Tianjin, Shijiazhuang, Tangshan, and Baoding, and had an outward expansion trend. For Chengde, Zhangjiakou in the north and northwest, and Cangzhou in the southeast, this type of area was mainly distributed in the city center. Moreover, since we mainly focused on the built-up areas of the 13 cities, which belonged to high-density human activity areas, there was almost no reverse urbanization. Therefore, the areas of type II and III areas were very tiny, and could be ignored. They are beyond the focus of this article.
Category IV refers to the decrease of the vegetation coverage in the process of urbanization, considering that the area initially used for plant growth was occupied by construction land. According to the above classification, the grid data were statistically analyzed at the pixel scale ( Figure 9). Due to the influence of outliers, the original statistical results made it difficult to show the partition. We thus excluded the exceptional values of two variables in order to show the distribution of the various pixels more intuitively.  Among the four categories, the proportion of Category Ⅳ was the highest, with a total area of 23,054.42 km 2 , accounting for 72.27% of the built-up area of the BTH region; the proportion of Category Ⅰ was the second, accounting for 26.89%, and 8281.73 km 2 ; the proportions of Category II (0.34%) and Category III (0.50%) were extremely tiny. This shows that the urbanization level in most of the built-up areas of BTH increased in the past 18 years, but there was apparent spatial heterogeneity in the vegetation growth in the different regions. With regard to Category I, it has two distinct distribution forms. On the one hand, it is distributed in the core areas of metropolises (e.g., Beijing, Tianjin, and Shijiazhuang). On the other hand, Category I was also located in the suburbs or coastal areas of Zhangjiakou, Chengde, Tianjin, and Hengshui. Category Ⅳ was mainly distributed in suburban areas, such as Beijing, Tianjin, Shijiazhuang, Tangshan, and Baoding, and had an outward expansion trend. For Chengde, Zhangjiakou in the north and northwest, and Cangzhou in the southeast, this type of area was mainly distributed in the city center. Moreover, since we mainly focused on the built-up areas of the 13 cities, which belonged to high-density human activity areas, there was almost no reverse urbanization. Therefore, the areas of type II and III areas were very tiny, and could be ignored. They are beyond the focus of this article.

Vegetation Variation and the Impacts of Urbanization on Vegetation
This study utilized the NTL index to characterize the level of urbanization by estimating the spatiotemporal dynamics of the FVC and NTL in the BTH region, in order to analyze the contributions of urbanization to FVC changes. We found that the vegetation coverage in the BTH region had apparent spatial heterogeneity. Due to the sufficient water resources in the northern Luan River Basin and the wide area of dense woodland, the FVC value of Chengde was higher than in other cities. As for the northwest of the BTH region and the Bohai Bay beach area, the FVC was relatively low. To be more specific, the low value of the FVC in the built-up areas was mainly attributed to the fact that the original vegetation cover was converted to impervious surfaces. As the precipitation in the BTH region was lower from 1999 to 2002 [53], the FVC value was lower from 2001 to 2003. Furthermore, after 2003, the precipitation increased; thus, the vegetation coverage rose. Moreover, the drought in 2009 contributed to the sharp decline in NDVI in 2009 [54]. From 2001 to 2018, the vegetation growth was enhanced significantly, with a pattern of greening in the northwest and browning in the southeast, which is consistent with the existing research conclusions [55]. The vegetation growth in the Zhangjiakou Bashang Plateau, the Yanshan-Taihang Mountains, and the eastern region of Cangzhou was noticeable, where the vegetation was mainly woodland, grassland, and shrubland [45]. Meanwhile, the vegetation of the central areas of Beijing, Tianjin, and Shijiazhuang showed significant growth as well. As for the areas where the vegetation coverage showed a significant downward trend, they were mainly located in the peripheral areas of the metropolises (i.e., Beijing, Tianjin, and Shijiazhuang) where there used to be a large area of cultivated land and impervious surfaces.

Differences in Responses of Vegetation Trends
Along the Urbanization Gradient, we discussed the mechanism of the impact of urbanization on vegetation dynamics in various categories, with the aim of providing a scientific basis for ecological protection and sustainable urbanization planning in the capital economic circle of China. We divided the areas into four major categories according to the changing trends of the NTL and vegetation dynamics in the Beijing-Tianjin-Hebei built-up area. We found that Category I and Category IV account for most of the areas of the BTH region, which account for 72.27% and 26.89%, respectively. Category I indicates that the vegetation coverage increased in the process of urban development, which means that there is a positive relationship between vegetation dynamics trends and urbanization. In the BTH region, a large amount of Category I areas are located in the core areas of metropolises (e.g., Beijing, Tianjin, and Shijiazhuang), which are more urbanized than other cities of the same agglomerations. In the advanced stage of urbanization, the negative effect of urbanization on vegetation will be countered by the positive demand for urban green space [56,57]. This theory can be verified in some developed countries, such as the United States, Britain, and Japan [40]. For Category IV areas, urbanization demonstrated a negative impact on vegetation, and these areas were mainly distributed in the suburban areas of metropolises. Meanwhile, the core urban areas of cities like Chengde, Zhangjiakou, and Cangzhou also belonged to Category IV. To be more specific, the suburbs of metropolises and the core areas of developing cities are in their early stage of urbanization. The rapid expansion of the city's original non-urban areas into urban areas resulted in a reduction in vegetation coverage [58]. This pattern of vegetation dynamics is common in developing regions, such as Africa [39] and the Middle East region [38].
Therefore, the spatiotemporal relations between urbanization and vegetation coverage dynamics are diversified, and they relate to the level of urbanization not only at the global, continent, and national scale, but also at the urban agglomeration scale. There are generally highly urbanized cities and developing cities in a specific urban agglomeration, which always demonstrates an apparent urbanization gradient. For the increasingly coordinated development of cities in the urban agglomeration, the vegetation dynamics generally show the following pattern: for the core urban areas of metropolises, whilst they are in the urban aggregation stage, continuing urbanization is effective for vegetation restoration, while for the periphery of metropolises and the core areas of developing cities that are still in the urban expansion stage, urbanization often has a negative impact on vegetation.

Effects of Ecological Engineering on Vegetation
We found that the vegetation in the suburbs of the Bashang Plateau area in the northwest and the coastal area in the southeast showed a significant greening trend, which may be attributed to the ecological engineering advocated by their authorities. Nowadays, the proposal of the 'Ecology Priority and Green Development Strategy' makes urban development gradually pay attention to the coordinated development of the ecology and economy, in order to alleviate ecological deterioration, which urged the government to implement environmental protection policies to improve environmental quality effectively.
The implementation of ecological engineering in ecologically fragile areas is a necessary procedure for urbanization, considering that the early stage of Chinese urban development neglected the importance of ecological protection and caused a series of environmental problems [59]. These environmental problems include the severe desertification in the northwest region and worse air quality in the downwind areas (e.g., Beijing and Tianjin). Consequently, the government has emphasized ecological protection in ecologically fragile areas and invested substantial funds in various ecological projects to ease the contradiction between ecology and economic development. Since 2000, a series of significant projects [60,61] such as the 'Grain for Green Project', the 'Beijing-Tianjin Sand Source Control Program', and the 'Three-North Shelterbelt Project' implemented in the north and northwest of the BTH region played a significant role in promoting vegetation protection.
We overlaid the boundary of the afforestation programs, the ecological function protection areas, and the FVC variation trend result of the BTH region from 2001 to 2018, in order to intuitively reflect the effects of various ecological projects ( Figure 10 Overall, the forestry projects and various ecological functional protection areas in the BTH region have not only alleviated the prominent contradiction of ecological problems (e.g., air pollution, water quality decline) caused by urban violent human activities [61,62] but also directly promoted the local vegetation greening. The ecological afforestation project in the sand source area directly improved the vegetation coverage in the sand source areas as well as alleviating the windblown sand problem in the downwind city [63]. Thus, with the continuous progress of human civilization, rapid urban development is effective in improving the awareness of human ecological protection and providing more funds and technology. Urbanization, to a certain extent, can help alleviate the environmental problems in the fragile ecological area and improve the vegetation cover at a particular stage.

The Correlation of Urban Greening and FVC Variation
For Chengde, Zhangjiakou, and Cangzhou, which were still in their early stage of urbanization, the vegetation growth in their core areas all showed a browning trend. This means that, with the increasing density of the urban center buildings, most cities faced the destruction of natural vegetation in their early stage of urbanization [64]. On the other hand, in terms of certain metropolises, the vegetation showed an apparent greening trend in the core areas of cities (Figure 8c). This difference may be attributed to the difference in the stages of urbanization between the urbanized cities and developing cities, suburbs, and rural areas [65]. That is to say, the effects of urbanization on vegetation varied during the different stages of urbanization [66]. For relatively highly developed areas, because of the increasing requirements from urban residents for the quality of living environments, the government emphasized the incorporation of the 'green coverage of the built district' as an assessment indicator of the cities' development level, which prompted Beijing, Tianjin, Shijiazhuang, and other similar metropolises to pay more attention to the vegetation growth in the built-up areas. The results of this study provide evidence that the effect of 'urban greening' in these cities is becoming visible. project in the sand source area directly improved the vegetation coverage in the sand source areas as well as alleviating the windblown sand problem in the downwind city [63]. Thus, with the continuous progress of human civilization, rapid urban development is effective in improving the awareness of human ecological protection and providing more funds and technology. Urbanization, to a certain extent, can help alleviate the environmental problems in the fragile ecological area and improve the vegetation cover at a particular stage. In this study, we chose Beijing, Tianjin, and Shijiazhuang as typical cities in order to analyze the relationship between vegetation coverage and the greening degree of built-up areas. According to the FVC trend graph of the BTH region, the urban core areas of the metropolises showed a definite 'green' trend from 2001 to 2018. We utilized the remote sensing monitoring data of land cover with a resolution of 30 m to obtain the built-up area range of the above three cities, and measured the mean FVC value of the built-up areas in 2001 and 2018. It was found that the average FVC value of the built districts in the three cities increased remarkably in the last 18 years. By referring to the China City Construction Statistical Yearbook from 2002 to 2018, we conducted a linear regression analysis of the green coverage in the built districts of Beijing, Tianjin, and Shijiazhuang, as shown in Figure 11. The green coverage of their built-up areas showed an apparent upward trend, which inferred that 'urban greening' had a positive effect on the greening of the vegetation coverage in the urban core areas. That is, in the highly urbanized areas, the relevance of the FVC variation trend and the urbanization level is positive, which is in agreement with previous studies [65]. It is due to the fact that metropolises can more easily be in pioneer positions to promote vegetation growth, especially in implementing advanced ecological restoration projects and advancing urban landscape planning. Sustainability 2020, 12,8550 Furthermore, the enhancement of vegetation in megacities can be attributed to many other factors. Urban greening plays multifunctional roles, such as the mitigation of climate change, the decrease of the atmospheric CO 2 concentration (and thereafter the urban heat island effect), and the increase of the air quality via decreasing active nitrogen [12,[67][68][69]. However, with the population and land resources showing apparent spatial mismatch, the maintenance of the balance between socioeconomic development and ecological sustainability still faces enormous challenges [70]. Government departments should promote the growth of vegetation according to particular local conditions and natural resource endowments.

Limitations and Further Improvements
Although nighttime light data have been prevalent in the studies of urbanization [71], there are still certain limitations on their usage. Firstly, the DMSP/OLS data itself has limitations in its coarse resolution and saturation effects, which bring uncertainties to the results [72]. For example, in this study, the results from highly developed cities, Beijing and Tianjin, have large oversaturation areas in the NTL index. Although a series of data corrections can mitigate the impact of the original data on the research results, it still leads to some deviation of the correction results of other areas from the actual light intensity, which has, to some extent, affected the measurement results of the NTL change trend. Therefore, in the subsequent research, it is necessary to focus on the methods of ameliorating the data correction in order to detect the development trend of the urbanization accurately. Secondly, NTL data cannot characterize the diversity, complexity, and dynamics of human activities completely. More research needs to consider diverse measurement factors in order to characterize the level of urbanization comprehensively. Besides this, in our study, there would be a slight influence on accuracy, which results from our approach to setting the DN value threshold for the purpose of obtaining the range of the built district of the BTH region. However, it serves well to reflect the scope of the built-up area, even though we realize that there are potential improvements.
There are a series of natural factors-such as climate [53], terrain [73], and groundwater [74]-that can affect regional vegetation growth. However, in this study, we mainly discussed the response of FVC in the BTH region to urbanization, considering that our foci are fractional vegetation and urbanization in metropolitan areas, and their interrelationships. Still, we have put these ecological factors into our discussion in order to interpret our results. Therefore, studies on the sensitivity of the vegetation growth of the BTH urban agglomeration to climatic factors (e.g., precipitation, temperature), seasonal drought, and other relevant factors are utterly promising. More research in the use of NTL data related to ecological, environmental, and social contexts should be undertaken, especially in metropolises, in order to explore the evolution of the ecological environment in urban agglomerations.

Conclusions
This study utilized the MODIS vegetation index product, DMSP/OLS nighttime lighting data, and statistical data to evaluate the spatial characteristic of vegetation coverage, the trend of the vegetation coverage variation, and the response of vegetation to the urbanization in a typical urban agglomeration, i.e., the BTH region in China.
The primary conclusions of this study include the following. First, the vegetation coverage in the BTH region had apparent spatial heterogeneity. From 2001 to 2018, the vegetation growth in the BTH region showed an upward trend, with a pattern of greening in the northwest and browning in the southeast. Second, the vegetation coverage was negatively correlated with the urbanization level. In the past 18 years, the increase in the intensity of human activity was accompanied by the destruction of local vegetation. Finally, in most areas (72%) of the built district in BTH, the vegetation tended to brown during urban development, while 27% of the built district was developing with improved urbanization and vegetation greening. Ecological programs and urban greening projects have played an essential role in the growth of vegetation in the BTH urban agglomeration. The ecological engineering in the ecologically fragile areas prompted the vegetation in the suburbs to green significantly; meanwhile, urban greening had a remarkable effect in the core areas of urbanized cities such as Beijing, Tianjin, and Shijiazhuang.
Overall, the vegetation's response to urbanization in a specific urban agglomeration shows apparent differences along the urbanization gradient, with the following pattern: for the core urban areas of metropolises, whilst they are in the urban aggregation stage, continuing urbanization is effective for the vegetation restoration, while for the periphery of metropolises and the core areas of developing cities that are still in the urban expansion stage, vegetation decreases in the process of urbanization. For urban agglomerations, urban greening projects can promote the growth of vegetation in metropolises, and other developing cities that are in their early stage of urbanization or are located in ecological fragile regions can promote their vegetation restoration by practical ecological projects. In the future, and the authorities of urban agglomerations have to consider the local urbanization gradient, and should formulate sustainable strategies scientifically in order to ensure the harmonious development of the economy and ecology.

Conflicts of Interest:
The authors declare no conflict of interest.