Spatiotemporal Patterns of Urbanization in the Three Most Developed Urban Agglomerations in China Based on Continuous Nighttime Light Data (2000–2018)

Abstract

Urban agglomeration is an advanced spatial form of integrating cities, resulting from the global urbanization of recent decades. Understanding spatiotemporal patterns and evolution is of great importance for improving urban agglomeration management. This study used continuous time-series NTL data from 2000 to 2018 combined with land-use images to investigate the spatiotemporal patterns of urbanization in the three most developed urban agglomerations in China over the past two decades: the Beijing–Tianjin–Hebei urban agglomeration (BTH), the Yangtze River Delta urban agglomeration (YRD), and the Guangdong–Hong Kong–Macao Greater Bay Area (GBA). The NTL intensity indexes, dynamic thresholds, extracted urban areas, and landscape metrics were synthetically used to facilitate the analysis. This study found that the urbanization process in the study areas could be divided into three stages: rapid urbanization in core cities from 2000 to 2010, a fluctuating urbanization process in both core cities and surrounding cities from 2010 to 2015, and stable urbanization, mainly in surrounding cities with a medium size after 2015. Meanwhile, the urbanization level of GBA was higher than that of YRD and BTH. However, with the acceleration of urban development in YRD, the gap in the urbanization level between GBA and YRD narrowed significantly in the third stage. In addition, this study confirmed that the scattered, medium-sized cities in YRD and GBA were more developed than those in BTH. This study showed that continuous NTL data could be effectively applied to monitor the urbanization patterns of urban agglomerations.

1. Introduction

Cities, as the centers of capital, labor, and information, were home to more than 55% of the world’s population in 2018, and this number is expected to reach 68% in 2050 [1]. Measuring spatial urbanization patterns is a crucial task in urbanization studies [2,3]. Angel and Blei [4] proposed the constrained dispersal structure of American cities based on the geographic information of workplaces and commuters. Taubenböck et al. [3] calculated spatial dispersion indexes and concluded that long-term urbanization created compact urban patterns. Under the processes of economic globalization and global urbanization, central cities and their neighbors develop to a specific spatial size with a high population density and continuity in the spatial landscape [5]. Cities, which have strong functional interconnectivity and accessibility, are highly integrated with a compact spatial organization and have gradually led to urban agglomerations [6]. These regions, which result from economic concentration and highly developed urbanized areas, have become the most important geographical unit for economic activities and technological innovations [7,8,9,10]. Compared with individual cities, urbanization processes in the urban agglomerations are more significant and intensified [11,12,13]. Hence, identifying the spatiotemporal patterns and characteristics of urbanization in such agglomerations is important for supporting sustainable development and management.

Based on the fact that cities are brightly lit at night, urban agglomerations can be easily distinguished in nighttime light (NTL) remote sensing images [14]. NTL data from the Defense Meteorological Satellite Program/Operational Linescan System (DMSP/OLS) have provided consistent global time series of night lights since the system’s launch in 1992 [15,16,17]. In 2013, the Visible Infrared Imaging Radiometer Suite sensor, fitted with a specific panchromatic sensor known as the Day/Night Band (VIIRS/DNB), succeeded DMSP/OLS as the monitor of global NTL [18]. Compared with annual land-use data or statistical data, NTL images can provide timely and precise information on the spatial extent of urban areas [14]. In addition, NTL intensity can indicate levels of urbanization, such as economic development and human activity intensity [19]. Hence, these two series of NTL datasets have been intensively used to retrieve data on the spatial and temporal dynamics of urbanization worldwide [20,21,22,23,24,25]. By utilizing multitemporal DMSP/OLS and VIIRS/DNB NTL annual composites, urbanization dynamics were observed in China [26,27,28,29], India [30], Japan [31], South America [32], and at global scales [33,34]. In recent papers, Georg et al. [35] spatially parameterized the characteristics of urban corridors and published a global inventory based on NTL images. Zhou et al. [36] developed a quantile-based approach to map the consistent global urban areas, finding that urban areas reached 0.53% of the world’s land surface area in 2013. Based on DMSP/OLS NTL data, Chen et al. [37] established a bidirectional Markov random field model and support vector machine classifier to extract urban areas and analyzed global urbanization progress from 2000 to 2012. Emerging studies on the development of urban agglomerations using NTL data have been reported in recent years. For example, Zou et al. [38] mapped the expansion in the middle reaches of the Yangtze River in China from 1992 to 2011 using DMSP/OLS data and concluded that, at different spatial scales and across temporal stages, the expansion directions and forms showed fluctuating characteristics of urban extents. Gao et al. [39] used the DMSP/OLS NTL dataset to validate the evolution of three urban agglomerations in China from 1992 to 2015. The result indicated that the city size distribution in the regions became more even, and the ranks of small cities ascended. Peng et al. [6] conducted a study that identified the boundaries of the 23 main urban agglomerations in China from 2000 to 2012 using NTL data. The study classified them as world level, regional level, or sub-regional level, which suggested the effectiveness of the NTL data in exploring the spatiotemporal process of urbanization in urban agglomerations.

Despite the widespread use of NTL data for monitoring spatiotemporal characteristics of the changing world [40,41], the lack of long-term time-series NTL observations limits investigations [42,43]. The existing studies generally focused on the urbanization development from 1992 to 2013 or from 2013 to the present [26,32,36], which was a consequence of the inconsistent temporal coverage of the two products. Although VIIRS/DNB products provide a higher spatial and radiometric resolution [44,45], as well as a wider dynamic range of anthropogenic light emissions [46,47], compared with previous products, the temporal duration is short, which limits its application in decadal observations. In addition, previously published studies mainly used NTL data to identify urban areas, impervious surfaces, and built-up areas [48,49,50,51,52], ignoring the distribution of light intensity and the internal characteristics of urban agglomerations [41].

The Beijing–Tianjin–Hebei (BTH), Yangtze River Delta (YRD), and the Guangdong–Hong Kong–Macao Greater Bay Area (GBA) urban agglomerations are China’s three largest urban agglomerations [39,53] and are the most populated regions in the world [54,55,56]. These three urban agglomerations have all experienced rapid urbanization, a large-scale expansion in urban population, and growth in NTL intensity [6,39,57]. The Chinese government has identified these three regions as large, national-level urban agglomerations, and the most important foci of China’s New-Type Urbanization Plan, which aims to accelerate their development and improve the coordination of city functions within them [5]. Strong economic development in these regions further encourages the government to engage in research on the spatial and temporal patterns of regional rapid urbanization. In general, BTH, YRD, and GBA clearly exemplify the country’s urban agglomerations that have experienced rapid urbanization.

Understanding the evolution of the spatiotemporal urbanization patterns not only provides policy guidance for these three regions but also provides a theoretical basis for urban management and improvements in regional competitiveness for other emerging urban agglomerations where urbanization is still in progress. However, few studies have focused on the spatial patterns of urban agglomeration development in China using long-term NTL images. Hence, based on extended time-series NTL data from 2000 to 2018, the spatiotemporal patterns of rapid urbanization in these three urban agglomerations were analyzed such that policymakers can propose sustainable and targeted strategies.

The objectives of this study were (1) to prolong the temporal range of urban agglomeration studies in China by using a continuous NTL dataset from 2000 to 2018 and (2) to identify the spatial dynamics and temporal evolution of the three most developed urban agglomerations. First, we calculated four NTL-related indexes to analyze the relationship between the NTL data and the level of urbanization, and to further find the metric that is most relevant for determining the internal dynamics of urban development. Second, optimal thresholds were determined for each NTL image to separate urban from non-urban areas for every year and to find threshold variation trends. Finally, a spatial, metrics-based analysis of urban areas that were extracted from the NTL data was conducted to reveal the overall spatiotemporal patterns.

2. Materials and Methods

2.1. Study Area

In our comparative analysis of urbanization spatiotemporal patterns, the three most populous and developed urban agglomerations (i.e., BTH, YRD, and GBA) were selected, as illustrated in Figure 1. These three agglomerations are also known as administrative urban areas, which were proposed by the Chinese government [58]. The total areas of BTH, YRD, and GBA are 218,000 km², 358,000 km², and 55,900 km², respectively. Their core cities are listed in

Table 1. Comparison of the socioeconomic statuses of the three most developed urban agglomerations in China in 2000 and 2018.

Urban Agglomeration	Total Area (km2)	Core Cities	Total Population (Million)	Percentage of Urban Population (%)	GDP (Billion USD)	GDP per Capita (USD)
The Beijing–Tianjin–Hebei urban agglomeration (BTH)	218,000	Beijing, Tianjin	107.57 (90.26)	65.84 (38.97)	1282.67 (124.06)	11,089.24 (1364.07)
The Yangtze River Delta urban agglomeration (YRD)	359,000	Shanghai	225.36 (195.53)	67.26 (42.51)	3195.81 (271.66)	14,229.61 (1375.56)
The Guangdong–Hong Kong–Macao Greater Bay Area (GBA)	56,000	Guangzhou, Shenzhen, Hong Kong	71.09 (49.99)	84.33 (72.90)	1653.96 (307.91)	24,310.08 (5908.35)

Note: The statistical data are included for both 2018 and 2000 (shown in the brackets). This statistical data can be obtained from the National Bureau of Statistics: http://www.stats.gov.cn/english/ (accessed on 7 June 2021).

. They have some similarities, such as a highly developed urbanization level, a megacity as the economic center, large populations, and a large GDP (Table 1). From 2000 to 2018, the residential populations increased from 90.26, 195.53, and 49.99 million to 107.57, 225.36, and 71.09 million, with GDP from 124.06, 271.66, and 307.91 billion USD to 1282.67, 3195.81, and 1653.96 billion USD, accounting for 9.23%, 23.00%, and 11.90% of the total GDP of China at the end of 2018 for BTH, YRD, and GBA, respectively. However, the development targets of these three agglomerations are significantly different due to their different population densities, political orientations, and economic conditions. In the future, BTH aims to be the most innovative region in China, with YRD and GBA as the pioneers of China’s economic reform and with an opening-up policy that aims to make them the most competitive urban agglomerations at the global scale.

2.2. Data Collection

This study used the extended time series (2000–2018) of the annual NTL images produced by Chen et al. [59]. This dataset, which presented good accuracy and temporal consistency, was generated using a modified auto-encoder model and a cross-sensor calibration model with two datasets: DMSP/OLS data (2000–2012) and composited monthly NPP/VIIRS data (2013–2018). The data, with a spatial resolution of 15 arc-seconds (approximately 500 m near the equator), using the unit of nW·cm⁻²·sr⁻¹, include stable night light from urban areas and towns. The instantaneous light caused by transient events, such as forest fires, was filtered. In addition, global land-use/cover data from 2000 to 2018, which were provided by the European Space Agency (ESA) Climate Change Initiative Land Cover (CCI-LC) dataset with a finer resolution of 300 m, was further used to identify local optimized thresholds as ancillary data in the study [60] due to its high overall accuracy among seven datasets of global land cover China [61]. All the remote sensing data were preprocessed to the Albers equal-area conic projection.

2.3. Methodology

The technical approaches used in this study can be broken down into three steps. First, an optimal NTL index was constructed based on the extended time-series NTL images to depict the urbanization level in each year of each urban agglomeration. Second, optimal NTL thresholds were determined using the ancillary data to accurately extract the urban areas. Finally, the evolution patterns of urban development were analyzed based on the NTL intensity.

2.3.1. Construction of an Optimal Nighttime Light Index

The variations in NTL intensity and lit areas can help to identify important characteristics of different urbanization stages. NTL intensity is a widely used indicator to estimate urban development, and a larger NTL intensity usually represents a greater level of development [62]. Lit areas, which are used for extracting urban areas, are well documented in mapping urbanization processes [14]. Compared with only using lit-area metrics, the index combining both NTL intensity and area, on the one hand, can reflect more detailed spatiotemporal patterns of rapid urbanization. On the other hand, an optimal index from remote sensing images can overcome the lag of statistical data and can be quickly obtained at the microscale. Therefore, the construction of an NTL index is an important link in the study of urbanization processes. To find an optimal index, this study first calculated four widely used indexes, of which the effectiveness has been shown in monitoring urbanization processes and urban expansion in China [58,63,64,65]. The correlations between these four NTL indexes and the economy, population, and land use of the three urban agglomerations were analyzed to determine the optimal index for analyzing the spatiotemporal patterns of the study areas. The first index was the luminous area, which is defined as the pixels with an NTL intensity larger than 0 nW·cm⁻²·sr⁻¹. The second index was the average NTL intensity index NTL_AI:

$${NTL}_{AI}={\displaystyle \sum _{r=\min }^{\max }r}\times \frac{nr}{N}$$

(1)

where r denotes the radiance of the pixels, n_r is the number of pixels with radiance r, and N is the total number of pixels in one image. This index took both light areas and intensity into account to indicate the urban internal characteristics.

The third index was the total NTL index:

$${NTL}_{TI}={\displaystyle \sum _{r=\min }^{\max }r}\times nr$$

(2)

where NTL_TI is the total NTL intensity of one image, r denotes the radiance of the pixels, and n_r is the number of pixels with radiance r.

The fourth index was a composited NTL index (CNLI), which was proposed by Zhuo et al. [65]:

$$CNLI={\displaystyle \sum _{r=\min }^{\max }\frac{r}{R_{\max }}}\times \frac{nr}{N}$$

(3)

where r is the radiance of the pixels, R_max is the maximum radiance of the image, nr is the number of pixels with radiance r, and N is the total number of pixels in one image.

After obtaining these four indexes, a Pearson correlation matrix was then established to examine the strengths of the bivariate associations between the indexes and statistical data, including economic levels, populations, and land uses.

2.3.2. Extracting Urban Areas by Determining Dynamic Thresholds of the NTL Image

Previous studies have proposed many approaches to separate urban areas from their background, such as Markov random field model [37], quantile-based approach [36], fully convolutional network [66], and support-vector-machine-based region [67]. This study used dynamic thresholds to extract urban areas for two reasons. First, dynamic thresholds do not ignore small cities or overestimate actual urban areas. Second, the variations in these thresholds indicate spatiotemporal changes in the urbanization level of the study areas.

For determining dynamic thresholds, this study integrated NTL images with CCI-LC land-use data. By comparing the total area of pixels with a radiance larger than i-class and the area extracted from land-use data, optimal thresholds for each NTL image were found when the difference reached a minimum value, as illustrated in Figure 2, where NI_max and NI_min are the maximum and minimum values of the NTL intensity, and their initial values are derived from the NTL image, NI_T is the threshold for the specific NTL image, NI_T-1 represents the value of NTL intensity that is the closest but smaller than NI_T, and NI_T+1 represents the value of NTL intensity that is the closest but larger than NI_T, S(NI_T) denotes the urban areas extracted from the NTL image using the threshold NI_T, Area is the actual urban area derived from the CCI-LC land-use data, and ΔS represents the discrepancies between urban areas extracted using threshold NI_T and Area values. The equation indicates that only when the threshold is NI_T can the discrepancy reach the minimum value.

After determining the thresholds, the urban area can be extracted. Considering that the land-use change within a single year was slight, we mapped the urban area with a five-year temporal interval, which is appropriate for monitoring urban expansion [68].

2.3.3. Evolution Patterns of Urbanization

The standard deviation ellipse (SDE) is an effective method for measuring the general trend and pattern of a changing phenomenon [69]. A change in the SDE can represent the evolution patterns of a discrete urban area, such as the change in gravity center, distribution range, density, direction, and shape, over the last two decades. The average center of the image extracted in Section 2.3.2 was a starting point to calculate one standard deviation of the x- and y-coordinates, and the long axis and the short axis of the ellipse were then defined [56]. The SDE was calculated using

$${SDE}_X=\sqrt{\frac{{\displaystyle {\sum }_{i=1}^n{(x_i-\overline{X})}^2}}{\mathrm{n}}}\mathrm{and}{SDE}_y=\sqrt{\frac{{\displaystyle {\sum }_{i=1}^n{(y_i-\overline{Y})}^2}}{\mathrm{n}}}$$

(4)

The center of the ellipse represents the relative position of the barycenter for the spatial distribution of urban areas, the long and short axes denote the dispersion degree of the two directions of the element, and the area of the ellipse represents the concentration degree of the spatial distribution of urban areas.

The intensity and compactness of urbanization development were measured using a relative expansion rate and cohesion metric, respectively. The relative expansion rate refers to the standardization of the average annual expansion rate of the urban areas in the three urban agglomerations, where this has been used to effectively quantify urbanization development [70,71]. The cohesion index is a measure of the compactness of an urban area [72,73], which increases as a patch becomes more clumped or aggregated in its spatial distribution.

$$R_i=\frac{A_i-A_{i-1}}{A_{i-1}}\times 100\%$$

(5)

$$C=\left(1-\frac{{\displaystyle \sum _{j=1}^n{p}_{ij}}}{{\displaystyle \sum _{j=1}^n{p}_{ij}}\times \sqrt{a_{ij}}}\right)\times \left(1-\frac{1}{\sqrt{N}}\right)\times 100$$

(6)

where R_i is the relative expansion rate of the urban area, i is the year, j is the serial number of the patch, A_i and A_i-1 are the area values of the urban patches in year i and i-1, p_ij is the perimeter of patch ij, a_ij is the area of patch ij, and N is the number of patches in the landscape. This class-level index was calculated using FRAGSTATS v4.2.

3. Results

3.1. The Change in NTL Intensity during Rapid Urbanization

Before analyzing the spatiotemporal change in NTL intensity, this study constructed a Pearson correlation matrix to determine the optimal NTL index. The Pearson correlation coefficients (

Table 2. Pearson correlation coefficients between the four NTL indexes and the urban statistical data of 2018.

NTL Index	Pearson Correlation Coefficients (Significance: Two-Tailed)
	Economic Metrics			Population Metrics			Land-Use Metrics
	GDP	GDP per Capita	Proportion of Non-Agricultural Industry Product in GDP	Population Density	Urban Population Density	Ratio of Urban Population	Construction Land Area	Built-Up Area
Luminous Area	0.067	0.268 *	0.088	0.725	0.329	0.068	0.507	0.611
	(0.599)	(0.028)	(0.634)	(0.096)	(0.635)	(0.836)	(0.432)	(0.371)
Average NTL intensity	0.928 **	0.896 **	0.593 **	0.851 **	0.643 **	0.379 *	0.753 *	0.803 **
	(0.000)	(0.001)	(0.000)	(0.000)	(0.000)	(0.042)	(0.017)	(0.003)
Total NTL intensity	0.830 **	0.892 *	0.482 *	0.809 **	0.631 **	0.436	0.648 *	0.705 **
	(0.000)	(0.028)	(0.023)	(0.006)	(0.000)	(0.129)	(0.043)	(0.000)
Composited NTL index	0.855 **	0.830 **	0.477	0.793 **	0.603	0.357	0.681 *	0.763
	(0.000)	(0.000)	(0.063)	(0.000)	(0.227)	(0.085)	(0.048)	(0.067)

Note: The result included both the coefficients (shown in bold) and significance (included in the parentheses). The statistical data can be obtained from the National Bureau of Statistics: http://www.stats.gov.cn/english/ (accessed on 7 June 2021). * The correlation was significant at the 0.05 level (two-tailed). ** The correlation was significant at the 0.01 level (two-tailed).

) showed that the four NTL indexes, except for the luminous area, were significantly related to the urban socio-economical characteristics. By comparing the coefficients of the NTL indexes, this study found that the average NTL intensity (NTL_AI) had the largest coefficients, which indicated that NTL_AI was the most optimal index for indicating the urbanization level and spatiotemporal patterns of urban agglomerations. Hence, this study selected NTL_AI as the light index for further analyses.

From 2000 to 2018, the NTL_AI continuously increased in the three urban agglomerations, which corresponded with rapid urbanization and economic growth. However, a significant variation of NLTAI in GBA and other two urban agglomerations was observed in this study (Figure 3D). This phenomenon could be partially attributed to the differences in the regional area to a much broader extent in BTH and YRD compared with GBA. Meanwhile, from the shape of the waterfall lines in Figure 3A–C, it is clear that the relative count of bright pixels in GBA was larger than that in BTH and YRD, which also caused obvious discrepancies in NTL_AI. After 2015, although the overall intensity maintained a relatively steady increasing trend, the pixels with low intensity in the three urban agglomerations experienced dramatic growth, which was attributed to development in rural and suburban areas.

3.2. Variations in the Thresholds for Extracting Urban Areas

This study calculated the yearly thresholds for the three urban agglomerations to accurately delineate the urban areas. The values of the thresholds can directly indicate urbanization level and socioeconomic development status across regions [74]. Generally, a higher threshold of NTL data represents a more developed urbanization stage [75]. The thresholds of the three urban agglomerations increased from 8.94, 2.95, and 1.92 in 2000 to 11.54, 9.95, and 2.92 in 2018 for GBA, YRD, and BTH, respectively. The continuous growth in the thresholds showed that the urbanization level in the three urban agglomerations kept increasing over the last two decades, which enhanced the urban area’s minimum luminous grade derived from higher-level living conditions. The result also indicated that the thresholds of GBA from 2000 to 2018 were the largest, and BTH ranked last (Figure 4). The highest thresholds and the NTL_AI of GBA suggested that GBA’s development level was much higher than the other two regions, revealing obvious discrepancies in the urbanization processes between the three agglomerations. However, it is worth noting that the slope of the trendline for the thresholds of YRD was much larger (0.32), i.e., three times and nine times larger than that of GBA (0.11) and BTH (0.037), respectively. Although the thresholds in YRD were only about 3 nW·cm⁻²·sr⁻¹ at the beginning of the study period, the acceleration of urban development made them very close to those of GBA in 2018. By contrast, the slope of BTH was relatively small with a low intercept, indicating that the urbanization development of BTH seemed to fall behind the other two regions [6,11].

3.3. Urban Areas Derived from NTL Images

Considering that the differences between yearly images were very small, this study selected NTL images for every five years (Figure 5). The NTL intensity distribution indicated that the three urban agglomerations had considerably high expansion magnitudes over the past two decades. It was apparent in Figure 5 that the pixels of core cities—such as Beijing and Tianjin in BTH; Shanghai in YRD; Guangzhou, Shenzhen, and Hong Kong in GBA—had high brightness. These core cities were so-called megacities and the economic centers of each urban agglomeration such that they played the most important roles in promoting the whole region’s development. Examining the distribution of high-level lit pixels indicated that the megacities of BTH expanded outward from urban to rural areas with ring structures, but the luminous area in YRD occurred both in the megacity, namely, Shanghai, and some subordinate cities, mostly in the Jiangsu and Zhejiang provinces. Megacities in GBA stood out with sector patterns due to terrain constraints. The distributions of lit pixels with high brightness in the three agglomerations were also observed in previous studies, which mainly focused on impervious surface expansion [53]. In general, the spatiotemporal patterns of the NTL intensity in BTH presented a rapid diffusion from the core areas of the central cities to the periphery with a pie shape, and the high-level lit pixels of BTH appeared much more compact than those in the other two regions. Meanwhile, the patterns in YRD and GBA were dominated by the development of scattered satellite cities in the outskirts. The differences might result from the particularly strong attractiveness of megacities in BTH, and small cities in the other two regions were more developed than those of BTH.

After determining the dynamic thresholds, this study extracted the urban areas of the three urban agglomerations on the basis of extended time-series NTL data (Figure 6). As mentioned before, the development pattern of BTH was driven by two megacities, with quite a few small cities. However, the urban area in GBA was much more scattered and presented a widespread trend from north to south [53]. The relative expansion rate of the urban areas suggested that the urban expansion in BTH and GBA experienced a similar trend. The trends both started at a high speed with a gradually decreasing expansion rate from 2000 to approximately 2010; during this period the rate decreased from over 12% to less than 5% in the two urban agglomerations. The curve then remained stable without any significant change. However, a difference occurred in the year 2015: the trend of GBA returned to an increasing state, while that of BTH remained at a low level. In contrast, the relative expansion rate in YRD showed obvious differences. Although the initial expansion rate was not very high, with a value less than 6%, the rate consistently increased from 2000 to 2013 and reached a peak of 14% in 2013. The expansion rate then rapidly decreased to about 3% in 2015 and the trend showed a slight increase after 2015. From the perspective of urban extents among the three urban agglomerations, GBA was the smallest but most developed region. Its urban area increased from approximately 4000 km² in 2000 to 10,000 km² in 2018, with an annual average growth rate of 2.4%. YRD, as the fastest-growing urban agglomerations in China, had the largest increase from 8500 km² in 2000 to more than 35,000 km² in 2018, and the annual growth rate was up to 6.5%. The expansion of BTH, which was mostly driven by two megacities, made the urban area increase at a speed of 3.9% annually from 8800 km² in 2000 to 26,000 km² in 2018. The urban expansion patterns indicate that the urbanization process in GBA was undergoing dramatic landscape changes since 2000, but the magnitude of the change was still much slower than that in the other two urban agglomerations.

4. Discussion

4.1. Evolution Patterns of Urban Areas

To further measure the evolution patterns and spatial directions of urban areas in the three agglomerations studied, the standard deviation ellipse and cohesion index of the urban expansion were calculated (Figure 7). The direction and concentration degree of the extracted urban area change at the regional scale were depicted by the central coordinate and area the ellipse. The spatial evolution of the three agglomerations presented significantly different patterns in northeast–southwest, northwest–southeast, and east–west directions in BTH, YRD, and GBA, respectively. The different directions of the long axis suggest the varied cluster modes of the three urban agglomerations. Considering the close distance between Beijing and Tianjin, the two focal points were exactly at two regional centers, with one at the capital of Hebei province and the other at the midpoint of the Beijing–Tianjin area. While satellite cities in YRD and GBA, such as Nanjing (YRD), Hangzhou (YRD), and Foshan (GBA), had strong economic vitality and intensive human activities, the direction of the long axis in these two agglomerations developed along the belt of these subordinate cities. Regarding the expansion or contraction of the two axes in the three agglomerations, the short axis expanded to a greater extent than did the long axis in BTH, indicating that the main driving force for the evolution of the urban area was the growth of cities along the northwest–southeast direction (the Beijing–Tianjin area) rather than the northeast–southwest direction. In contrast, the rapid urbanization of peripheral cities in GBA made the long axis of the ellipse expand to a greater extent than the short axis. However, the two axes in YRD both contracted during the period, which suggests that the key drivers were cities located around the ellipse’s center. These findings corroborate the evidence reported in previous studies [58,76], confirming that scattered medium-sized cities in YRD and GBA were more developed than those in BTH.

Table 3. The parameters from the directional distribution analysis using standard deviation ellipses and cohesion indexes of urban areas in three urban agglomerations for 2000, 2005, 2010, 2015, and 2018.

Region	Year	Short Axis (X) (km)	Long Axis (Y) (km)	Latitude of Center Point	Longitude of Center Point	Rotation	Cohesion Index
BTH	2000	97.95	219.23	39°03′24″	116°15′21″	23.53	90.12
	2005	98.02	217.88	38°59′06″	116°14′50″	23.00	90.80
	2010	100.06	218.47	38°55′15″	116°16′04″	21.83	90.86
	2015	101.96	220.62	38°53′00″	116°16′53″	20.72	89.91
	2018	103.67	220.26	38°53′12″	116°18′30″	20.98	91.25
YRD	2000	151.17	317.22	31°15′10″	119°48′34″	139.43	87.17
	2005	154.33	330.65	31°29′16″	119°30′54″	139.10	86.24
	2010	151.84	320.33	31°33′23″	119°31′57″	138.80	89.18
	2015	146.61	314.32	31°32′07″	119°38′08″	138.57	91.62
	2018	147.85	309.61	31°31′41″	119°37′44″	138.35	92.83
GBA	2000	72.72	45.65	22°51′45″	113°31′40″	96.82	93.88
	2005	71.34	46.50	22°52′17″	113°34′32″	92.66	97.37
	2010	71.33	47.90	22°52′41″	113°32′04″	91.74	98.52
	2015	71.77	49.11	22°52′47″	113°30′51″	91.21	98.94
	2018	71.59	49.74	22°52′44″	113°30′11″	89.70	99.24

shows that the rotation of the ellipse of BTH, YRD, and GBA decreased from 23.53, 139.43, and 96.82 in 2000 to 20.98, 138.35, and 89.70 in 2018, respectively, which meant that the spatial patterns of the urban areas gradually changed from the original direction to new statuses with the rapid development of urbanization. Based on the changes in the cohesion indexes, the compactness of the three regions seemed to increase to different extents. Among them, the cohesion of BTH experienced the slowest increase over the past two decades, with the smallest value in 2018. YRD and GBA, however, both had significant growth and reached a high level of connectivity. This discrepancy could be attributed to how closely the inner cities interacted with each other. Small- and medium-sized cities in BTH lacked economic communication, limiting the cohesion of the whole region.

In fact, these significant differences in the evolution of the urban areas in the three urban agglomerations were one of the principal manifestations of spatiotemporal variations in the rapid urbanization progress. The urban areas extracted from the long-term NTL images consistently provided us with deeper insight into these variations at a low cost.

4.2. Analysis of Spatiotemporal Patterns of Urbanization Based on Multiple Metrics

The NTL_AI, dynamic thresholds, SDE, cohesion index, and characteristics of urban area expansion played important roles in understanding the spatiotemporal patterns of the rapid urbanization process in the three urban agglomerations. Combined with the multiple metrics, it can be concluded that these urbanization processes could be divided into three stages.

The first stage was the rapid development of urbanization, with consistently increasing human activity intensity and significant urban central area expansion from 2000 to 2010. During this period, the relative expansion rate was larger than 3%, with NTL_AI and the thresholds presenting obvious increasing trends. Meanwhile, the migration of the center points and the decrease in the rotation for the ellipses were significant (Figure 7 and Table 3). Considering that the outline of China’s Tenth and Eleventh Five-Year Plans emphasize the development of core cities in urban agglomerations, a possible explanation for these phenomena is that the rapid expansion mainly occurred in the core cities rather than the surrounding municipalities. Strong growth in both land-use intensity and economic activities, with population migration in the core areas of megacities, made the urban agglomeration enter its rapid urbanization stage. At this stage, the main manifestation of NTL was the drastic intensity increase in local areas such as Beijing, Shanghai, Shenzhen, and Guangzhou. Due to the increase in the overall NTL intensity, dynamic thresholds experienced significant growth. The rapid urbanization of megacities intensively promoted the whole region’s development. Previous studies showed that coastal areas exhibited the most dramatic growth rate in urban areas [68] and a higher urbanization level [6] in this period. The urbanization of GBA was much more developed than the other two regions, with a higher GDP per capita, more intensified NTL, larger thresholds, and a higher cohesion index.

In the second stage, the urbanization level increased at a fluctuating rate, with obvious discrepancies between the three urban agglomerations, from 2010 to 2015. Since the urban development in GBA was ahead of the other two regions, the urbanization progress in GBA significantly slowed down at this stage, while BTH and YRD were still in relatively high-speed development. Combining the variation in the NTL_AI and lit pixel count with spatial patterns of urban expansion, this study found that human activities and urban expansion in YRD and GBA gradually shifted to satellite cities, such as Nanjing (YRD), Hangzhou (YRD), and Foshan (GBA), which might have resulted from the Strategies of Regional Coordinated Development proposed by the Chinese government. Due to the relative dimming of NTL in subordinate cities, the thresholds in this period fluctuated. However, the gravity center of development in BTH was still strongly affected by Beijing and Tianjin. At this stage, the sizes and magnitudes, as well as the spatiotemporal patterns of the urban areas, varied significantly. GBA and YRD gradually formed into megapolises with scattered urban development patterns and subordinate cities or satellite cities sprawled outward, while BTH grew mainly due to the impetus from core cities.

In the third stage, areas still undergoing the urbanization process mainly transformed from core cities to satellite cities in the three urban agglomerations with a stable expansion rate after 2015. China’s Thirteenth Five-Year Plan accelerated the construction of medium-sized cities in the urban agglomerations studied to optimize the regional structure. It was apparent that the count of lit pixels with low NTL intensity, compared with the high-brightness pixels, increased much more rapidly after 2015 (Figure 3). Meanwhile, the thresholds of NTL in the three agglomerations slightly increased. These signals indicated that urbanization in the three agglomerations focused on subordinated cities, even in BTH, which was relatively less developed. At this stage, quite a large number of medium-sized cities underwent urbanization processes due to the “saturated status” in major megacities [77]. Hence, the relatively slow urbanization development in these surrounding cities led to only a slight increase in the NTL_AI, thresholds, and expansion rate in the agglomerations. This study also found that NTL_AI in YRD and BTH was much lower than that of GBA in the three stages. A recent paper provided outlines of the urban agglomerations in a consistent, data-driven way and reported that the morphological urban area of GBA was the largest in the world, while that of YRD and BTH ranked at numbers 10 and 17, respectively [78]. Consequently, the significant differences in NTL_AI between the three agglomerations could presumably be attributed to a much larger discrepancy between administrative areas and morphological areas in BTH and YRD compared with GBA.

4.3. Limitations and Future Study

This research was conducted for BTH, YRD, and GBA, which are typical urban agglomerations in China. This study defined urban agglomerations using their administrative boundaries; however, Taubenböck et al. stated that the actual morphological urban areas are smaller than the administrative units [78]. This requires further consideration in future studies. Although previous studies pointed out that YRD met the world-level standard, discrepancies between the three urban agglomerations and the most developed regions of the entire globe were shown [53]. Comparison studies on the urbanization processes at the global scale using long-term time-series observation data are recommended. This study used extended time-series NTL data without any comparison with other products, such as Luojia satellite images [79,80] and ground-based measurements [81]. Hence, using composited NTL data from multiple sources might assure the robustness of consistent time-series images. In addition, China’s “ghost neighborhoods” phenomenon has been identified in dynamic regions, such as YRD and GBA [82], which might underestimate the NTL intensity. Therefore, further studies that combine NTL images with real-world activity intensity are desirable.

5. Conclusions

This study conducted a multi-metric-based analysis of the three most developed urban agglomerations in China from 2000 to 2018 based on long-term time-series NTL data by characterizing the spatiotemporal patterns of rapid urbanization in the regions. The three urban agglomerations, namely, BTH, YRD, and GBA, had significantly different urbanization patterns and rates. In general, the urbanization level of GBA was the highest among the study areas. Not only the high NTL intensity but also the large urban area and the developed satellite cities indicated that GBA was the most developed urban agglomeration in China. However, the economic development and rapid population growth dramatically and intensively accelerated the urbanization processes in YRD, of which the NTL_AI almost surpassed that of GBA. The satellite cities in both GBA and YRD entered a stage of rapid expansion and high connectivity with the megacities of each urban agglomeration. The development of BTH was still very competitive due to the strong impetus from Beijing and Tianjin, but peripheral cities should be paid more attention to regional urbanization.

Based on the spatiotemporal patterns of urban development in the study areas, the urbanization processes can be divided into three stages. In the first stage, from 2000 to 2010, all three urban agglomerations experienced rapid urbanization, including a significant increase in NTL_AI and thresholds, and a rapid expansion of the urban central areas. In the second stage, from 2010 to 2015, the urbanization level increased at a fluctuating rate, with emerging discrepancies between the three urban agglomerations. The peripheral cities in GBA and YRD played increasingly important roles in regional development, while the urbanization of BTH was still driven by megacities. The third stage, starting in 2015, was the developed stage. The spatial distribution of areas where urbanization was still in process mainly moved from megacities to medium-sized cities with a relatively low NTL intensity. A comparative analysis of spatiotemporal patterns of urbanization in the three urban agglomerations might provide valuable lessons for urban management in other regions at different development phases.