Abstract

Impervious surface (IS) is a key indicator to measure the urbanization process and ecological environment. Many studies have observed an urbanization process based on IS at the city scale. Understanding the changes in the IS over a period at a regional level offers an alternative and effective approach to characterize and quantify the spatial process of urban agglomeration. This study focuses on the urban agglomeration of the Guangdong-Hong Kong-Macao Greater Bay Area (GBA) by utilizing the sensor-based Landsat data during 1987-2017 and investigates the spatiotemporal distribution of IS expansion at both regional and city scales. The modified linear spectral mixture analysis (MLSMA) method is used to extract the IS of the GBA. Then, the IS mapping accuracies were assessed after comparison with high-resolution historical data. The spatiotemporal and directional changes of IS surfaces for GBA are analyzed by using Gravity Center (GC) and Standard Deviational Ellipse (SDE). Finally, Shannon’s Diversity Index (SHDI) is used to analyze the overall characteristics of landscape level, and the Patch Density (PD) and Landscape Shape Index (LSI) are used to describe the characteristics of different classes of the IS. The results show that the IS of the whole region experienced rapid and massive expansion during the past 30 years and exhibited a distinct characteristic along the Pearl River Delta (PRD) and the coastline. Furthermore, the IS area increased rapidly in the PRD, while it is relatively stable in Hong Kong and Macao. We believe that the findings of this study can help policy makers to better understand and maintain the sustainable development of the GBA.

1. Introduction

In 1950, 30% of the world population resided in the urban areas due to global urbanization which has progressed with unprecedented speed and thus resulted in the increased urban population to 54% in 2014 [1]. With this global urbanization, in the past two decades, China has also experienced rapid urbanization which replaced natural vegetation coverage with impervious surface (IS), through which water cannot infiltrate into the soil such as buildings, paved roads, and parking lots [25]. Therefore, IS change indicates urban development and can reveal significant information about urban areas; they can be utilized to quantify urban development and the spatial process of urban agglomeration. Way back in 1898, Howard [6] proposed a concept of “Town Cluster,” followed by global researchers and opened debate on concepts like “urban economic region,” “economy city,” conurbation,” and “megalopolis” [79]. In the lights of the aforementioned concepts, Chinese scholars put forward the concept of “urban agglomeration” [10, 11]. An urban agglomeration is defined as a type of megalopolis [12], which consists of one or more megacities within a core area with at least three large cities within a specific geographical area. The development of the core areas is measured with the advancement in the transportation system, communication services, integrated economic structure, and socio-culture functions. Thus, the urban agglomeration spatial distribution characteristic reflects the internal spatial structure [13] and the gaps in the economic development of the core areas. However, the urban evolution and ecological construction of different urban agglomerations are mainly caused by the differences in spatial structures [1418] of the cities. Therefore, studying the spatial structure of urban agglomeration with the combination of IS can be of great significance to promote spatial planning and development.

Since the 1980s, a substantial amount of research has focused on understanding the urban agglomeration spatial structure to cover issues such as the spatial expansion mode [19, 20], evolution character [21, 22], and driving mechanism [2325]. Lu et al. [26] examined the land expansion of urban agglomeration in whole and local regions of Wuhan, China, by using a series of urban expansion indexes and discussed its influencing factors as well as the problems. Gao et al. [27] revealed that the urban land of the Yangtze River Delta (YRD) experienced the growth of residential and industrial land and investigated the driving forces of urban land expansion and structural changes by multimodels. Man et al. [28] examined the characteristics of spatiotemporal IS variation, quantified the progression of the variation using IS change trajectories, and explored the orientation of IS expansion in the context of rapid urbanization. Ma et al. [29] investigated the impacts of IS area on land surface temperature (LST) in the urban agglomeration region and indicated that landscape metrics of IS and the density of IS had a significant correlation with the LST. Moreover, Yang et al. explored the spatiotemporal evolution of 11 cities within the Guangdong-Hong Kong-Macao Greater Bay Area (GBA) by integrating remote sensing, landscape analysis, and geographic information system (GIS) techniques and stated that the expansion of most cities took on an urban-to-rural landscape gradient [30]. Yang et al. [31] argued that IS could be used as a synthesized quantifiable index to reflect the intensity of natural ecosystems changing into urban ecosystems. Furthermore, the authors stated that the results improve the accuracy and efficiency of extracting IS and analyze the changing characteristics of IS in the context of China’s urbanization. Cao et al. [18] used a gravitational index for the comparative analysis of the spatial structure of IS between Chinese and American city clusters and indicated that most of the current researches were focused on the driving force and effect of the spatial structure evolution of urban agglomeration in whole and local regions. However, urban agglomeration contains different scales like urbans, metropolis, and metropolitan areas; the development and influence factors are much more complicated [3234].

Furthermore, understanding the long-term land expansion and spatial structure change in various geographical scales and locations is better to understand the historical development and the future trend. Therefore, it is necessary to investigate the spatiotemporal IS change of urban agglomeration in various scales. The IS greatly effects the regional climate [34, 35], environment [36], and the water resources [37] and is considered an important indicator for urban planning, environment, sustainability, and resource management [4, 5] and is utilized to the study of urban agglomeration spatial structure. Remote sensing is considered the primary method to extract IS with several advantages (i.e., multitemporality, broad coverage, and low cost) as compared to the traditional methods [38]. Recently, numerous studies have been conducted to investigate the methods of IS extraction, such as stepwise multivariate statistical method [39, 40], artificial neural network [41, 42], decision tree method [43, 44], spectral mixing model [45, 46], and spectrum-based index method [4751]. All the methods can effectively extract IS. However, the medium-resolution image is recommended and is widely adopted to monitor the process of urban spatial structure due to its continuous measurement at moderate spatial resolution and temporal frequency [52]. Currently, a combination of the V-I-S model and spectral mixture analysis based on the remote sensing data is considered the most popular method to estimate the urban IS [53, 54] as compared to the conventional [55]. In order to improve the classification results, a modified linear spectral mixture analysis (MLSMA) method was proposed by Xu et al. [56], which combines a Normalized Difference Built-up Index (NDBI) and Normalized Difference Vegetation Index (NDVI) with Linear Spectral Mixture Analysis Method (LSMM), and is widely used [57].

GBA is one of the major bay areas in the world, and it has been experiencing an intensive urbanization process since the reform and opening-up policy in 1978 [58, 59]. The motivation of choosing GBA as our research area is the cross-administrative distinction region, the scales of administration, and the complicated geography. We believe that the study on the IS changes for GBA not only provides a decision-making platform for the planning and management but can also provide a reference to better understand the spatial structure evolution of urban agglomeration. The purpose and innovation of this study are to introduce the methods of IS change on the urban agglomeration for GBA by utilizing the Landsat images based on landscape pattern index for 1987, 1992, 1997, 2002, 2007, 2012, and 2017 and to investigate the spatiotemporal distribution of IS expansion at both regional and city scales. In this study, the MLSMA method is used to extract the IS of the GBA by using sensor-based Landsat data (i.e., Thematic Mapper (TM) Sensor, Enhanced Thematic Mapper Plus (ETM+) Sensor, and Operational Land Imager (OLI) Sensor). Then, the IS mapping accuracies were assessed after comparison with high-resolution historical data. The spatiotemporal and directional changes of IS for GBA were analyzed by using Gravity Center (GC) and Standard Deviational Ellipse (SDE), which include the following steps: (a)Comparing the extension differences between PRD and two special administration regions by calculating the IS growth speed and rate(b)Quantifying the spatiotemporal structure evolution of three metropolitan areas in GBA by using the GC and SDE analysis

Finally, the similarities and differences of spatial and temporal change at the landscape pattern of core cities in GBA were compared by using Shannon’s Diversity Index (SHDI), Patch Density (PD), and Landscape Shape Index (LSI).

The organization of the rest of the paper is as follows: Section 2 defines the study area and dataset. Section 3 presents the methodologies used. Section 4 presents the results and discussion for the experimental results performed on the dataset. Finally, Section 5 concludes the study and proposes a future work.

2. Study Area and Dataset

2.1. Study Area

To strengthen the cooperation between the Chinese mainland and Hong Kong-Macao and establish a new development pattern, GBA was proposed in the government work report of the 12th National People’s Congress [60] by the urban planners and government. GBA is an urban agglomeration formed by the combination of nine cities (Guangzhou, Shenzhen, Zhuhai, Foshan, Huizhou, Dongguan, Zhongshan, Jiangmen, and Zhaoqing) in the Guangdong province and the two Special Administrative Regions Hong Kong and Macao, as shown in Figure 1. It forms the “” economic spatial structure “Guang-Fo-Zhao” (Guangzhou, Foshan, and Zhaoqing), “Shen-Guan-Hui” (Shenzhen, Dongguan, and Huizhou) and “Zhu-Zhong-Jiang” (Zhuhai, Zhongshan, and Jiangmen). Meanwhile, it forms three development poles, which include “Guangzhou-Foshan,” “Hong Kong-Shenzhen,” and “Macao-Zhuhai.” Guangzhou, Shenzhen, Hong Kong, and Macao are four core cities in this area. It is an important economic and industrial zone for China’s development of a world-class urban agglomeration and participation in global competition. It is also included in one of the four largest bay areas in the world (i.e., New York, San Francisco, and Tokyo). GBA covers 56000 km2, and the population grew to 66.72 million by the end of 2017, and the GDP of GBA reached 10 trillion RMB [61]. The economic aggregate of GBA has exceeded that of the San Francisco Bay Area and is almost equal to that of the New York Bay Area. As an urban agglomeration that acts as a carrier of China in global competition, the cities in GBA have experienced rapid urbanization; therefore, it has great significance to study spatiotemporal structure changes to better understand the sustainability and urban development of the GBA.

2.2. Dataset and Preprocessing

The remote sensing images utilized in the current study are acquired by Landsat 5 TM, Landsat 7 ETM+, and Landsat 8 OLI sensors over the period 1987 to 2017 with the 30 m spatial resolution. The administrative boundary of the GBA is located within eight scenes (path/row: 121/44, 121/45, 122/43, 122/44, 122/45, 123/43, 123/44, and 123/45). The datasets used in the study are the surface reflectance dataset provided by USGS [62]. All Landsat surface reflectance dataset for the study area is almost cloudless. The ground control points (GCPs) were geometrically correct for the remotely sensed images, and the square root error of the corrected image pixels is less than 0.5 pixels. The digital number (DN) value of the image brightness value is converted to the standard on-board reflectance by radiometric calibration to eliminate the difference in sunshine conditions in the multispectral image [63, 64]. Lu and Weng [63] concluded that the atmospheric correction does not contribute significantly to the extraction of the IS, so the remote sensing images utilized in this study were not corrected for the atmosphere.

In the current study, we used historical high-resolution images acquired from Google Earth [65] as reference data to assess the accuracy of the IS area. Besides, the topographic map of the Institute of Surveying and Mapping, Department of Natural Resources of Guangdong Province [66], the Digital Orthophoto Map (DOM) data of Guangzhou Jiantong Surveying, Mapping and Geoinfo Co., Ltd. [67], and field survey are used for verification. Meanwhile, the boundary file of administrative units (downloaded from the National Geomatics Center of China [68]) is used to analyze the spatial distribution of the IS area in the study area.

3. Methods

In the current study, a modified linear spectral mixture analysis (MLSMA) method proposed by Xu et al. [57] is used to extract the impervious surfaces of the GBA in 1987, 1992, 1997, 2002, 2007, 2012, and 2017. Based on the time-series impervious surface fraction images, the spatial distribution pattern and differentiation law of IS in the research area are analyzed by using spatial and landscape pattern analysis and revealed the process of surface cover change in the megalopolis development from 1987 to 2017.

3.1. IS Mapping
3.1.1. The Modified Linear Spectral Mixture Analysis Method

Mapping the dynamics of the impervious surface and measuring orientation and direction of IS in the study area is important to understand the spatial and temporal change characteristics of urban development during the past 30 years. The change of the IS area is mapped using the MLSMA in three steps, as shown in Figure 2. (a)Calculate the modified normalized difference water index (MNDWI) and then mask the water body of Landsat images by the result of MNDWI. The threshold values of the MNDWI for 56 Landsat images of different times were experimentally determined in this study(b)Extract the fractions of high and low albedo, vegetation, and soil by the conventional LSMM method(c)Extract the IS, vegetation, and soil fractions with the unmixing results of the conventional LSMM, Normalized Difference Built-up Index (NDBI), and Normalized Difference Vegetation Index (NDVI). The threshold values of NDBI for 56 Landsat images of different times were experimentally determined in this study. A previous study has shown that pixels with a NDVI value less than 0.2 are classified as bare soil [69]. Therefore, this study set pixels with NDVI values less than 0.2 as soil and greater than 0.2 as vegetation

3.1.2. Accuracy Assessment

The IS mapping accuracies were assessed after comparison with high-resolution historical data acquired from Google Earth, Institute of Surveying and Mapping, Department of Natural Resources of Guangdong Province and Guangzhou Jiantong Surveying, Mapping, and Geoinfo Co., Ltd. The available historical images from Google Earth are limited and cannot cover the whole study area before 2007. Thus, the data of 2012-2017 are used to assess the accuracy of the contemporary IS classification result. The three steps involved to assess the accuracy are as follows: Firstly, we carried out the space registration of reference and classified. Secondly, we collected 500 IS testing samples and 500 pervious surface testing samples in each city and the whole region randomly by stratified random sampling [64] rule. Lastly, we compared the results with visual interpretation and then calculated the overall accuracy, sensitivity, and specificity to evaluate the classified result [70, 71] by computing the confusion matrix [72] as follows: where the connotation and relationship of TP, TN, FP, and FN are shown in Figure 3.

The higher the value of sensitivity and specificity, the smaller the probability of misclassification and omission of IS, although the precision of IS during a few periods cannot be assessed due to the absence of high-resolution images. Thus, it is assumed that they are similar to the result of 2012-2017.

3.2. Spatial Directional Analysis

Detecting the spatiotemporal and directional changes of the IS can provide important information for optimizing the regional planning and management. In the current study, the spatiotemporal and directional changes of IS on three metropolitan areas Guang-Fo-Zhao, Shen-Guan-Hui, and Zhu-Zhong-Jiang, and three development poles Guangzhou-Foshan, Hong Kong-Shenzhen, and Macao-Zhuhai in the bay area are analyzed by using GC and SDE. The GC is employed to identify the weighted center of geographic elements [73], and the SDE is used to explain the centrality, distribution, directivity, and the spatial form characteristics of global space and reveals the temporal and spatial distribution of geographic elements accurately [74]. There are four parameters of GC and SDE, which include the gravity center, long and short axis, and azimuth.

The gravity center is calculated as follows:

The azimuth “” of SDE is calculated as follows:

The standard deviations of the ellipse “” and “” in directions “” and “” are calculated as follows: where is the spatial location of the object, is the corresponding weight, is the coordinates deviation from the spatial location, and () is the gravity center.

3.3. Landscape Pattern Analysis

By using the equidistant division [75], the IS fractions of the core cities (Guangzhou, Hong Kong, Macao, and Shenzhen) are segmented into five categories: very low density, low density, medium density, high density, and very high density, with the interval of (0%,20%], (20%,40%], (40%,60%], (60%,80%], and (80%,100%].

In the current study, the landscape pattern analysis is adopted to indicate the IS features and the difference in the biophysical composition of the same IS category to better understand the difference of spatiotemporal characteristics among the core cities. To analyze the overall characteristics of landscape level, SHDI is used, and the PD and LSI are used to describe the characteristics of different classes of the IS [78]. The implication of these indices is listed in Table 1.

4. Results and Discussion

4.1. Spatiotemporal Patterns of IS in Guangdong-Hong Kong-Macao Bay Area
4.1.1. The Overall Spatiotemporal Pattern Changes of IS in the Bay Area

Based on the sensor-based Landsat images, the distribution of the IS surface is extracted for 1987, 1992, 1997, 2002, 2007, 2012, and 2017. To calculate the overall accuracy, sensitivity, and specificity, high-resolution images from Google Earth are selected to assess the classified accuracy and the reliability of the IS mapping for 2012-2017. All accuracy metrics of the classification were greater than 85% in 2012 and 2017 (Table 2), which fulfilled the accuracy requirement. For verification due to the lack of historical high spatial resolution data for reference, data from topographic maps of Institute of Surveying and Mapping, Department of Natural Resources of Guangdong Province, Digital Orthophoto Map (DOM) data of Guangzhou Jiantong Surveying, Mapping, and Geoinfo Co., Ltd., and field survey are combined to meet the precision for analysis by considering the size of the study area.

Figure 4 presents the IS distribution in 1987, 1992, 1997, 2002, 2007, 2012, and 2017 and reviles that the IS increased significantly with the passage of time in the last 30 years, and it has the distribution characteristic which takes the Pearl River and coastline as the axis. In the development process of GBA, cities along the tributaries of the Pearl River are the core regions of development and construction. This phenomenon fits the universal law of urban construction and economic development and highlights the GBA’s policy of building economic belt along the river. Moreover, the intensity of urbanization in the coastal regions of GBA is greater than that in other regions at an early period, as these regions are the frontier regions of the reform and opening-up policy.

Figure 5 presents the IS area, growth rate, and speed of the IS of the bay area at different time periods. It can be observed from Figure 5(a) that from 1987 to 2017, the IS area increased from 1939.338 km2 to 12384.95 km2, with an annual average speed of 351.52 km2/a. From 2007 to 2012, the IS area experienced the fastest growth with a speed of 622.16 km2/a. However, Figures 5(b) and 5(c) reveal that the growth rate and speed of IS in GBA experienced “down-up-down” tendency, which indicates that the IS expansion of the GBA is under expansion with a relatively slow paced development since 1970’s rapid development [71] which has executed the rapid transition from agricultural landscape to metropolis [79].

4.1.2. The Spatiotemporal Pattern Changes of IS in the Cities of GBA

The IS area of GBA increased significantly (see Figures 6 and 7) from 1987 to 2017, while other cities in GBA have different temporal and spatial characteristics of IS expansion. Figure 6 presents the IS area, growth rate, and speed of the cities in GBA. Generally, the cities in GBA have expanded rapidly since 1987 except Hong Kong and Macao. The growth rate of IS in Hong Kong is relatively low due to the limitation of a geographical environment. The IS in Macao started to expand rapidly in the 1990s. As mainland cities in GBA started the rapid urban construction after the 21st century, IS area increased greatly. Among them, the development of Dongguan, Foshan, Guangzhou and Zhongshan was more remarkable. However, Shenzhen has experienced a rapid urbanization since 1987 as it is close to Hong Kong. The total areas of Guangzhou, Foshan, and Jiangmen are very large due to their plenty expand spaces and the strong driving force of regional economic development.

IS expansion of the cities in GBA is mapped in Figure 7. The result shows that the urban expansion of the cities in GBA has significant spatial and temporal differences.

In Hong Kong and Macao, the two special administration regions, as a geographical environment restricts the urban expansion, IS changing was dominated by internal adjustment. The development and expansion of IS mainly occurred from 1987 to 1997 and then developed steadily, and the structure of these two cities becomes stable. Moreover, the growth rate and speed of IS in these two cities are lower than those in the other nine cities in mainland China.

In Shenzhen and Zhuhai, the IS area was primarily distributed at the intersection with Hong Kong and Macao in the early period, then expanded to multicore. From 1987 to 1997, the IS increased with a high growth rate and speed, which in turn generated the core of the city. However, the growth rate and speed of IS began to decline after 2002, and the structure of Shenzhen tends to stabilize. Moreover, the IS increased at the intersection with Zhuhai during 1987-2002. From 2007 to 2017, the IS experienced internal growth and almost covered the whole city. As a special economic zone, the IS in Zhuhai mainly appeared at the intersection with Macao. The IS expanded significantly with a high growth rate and speed from 1987 to 2002.

Foshan and Guangzhou showed an obvious development trend of the core region driving the edge region. In Foshan, the growth rate and speed of the IS showed “M”-type development, while after 2002 it increased significantly as compared to the period of 1987-2002. Moreover, the IS was mainly distributed along the border between Guangzhou and Foshan in the early period, which then expanded to the whole city. From 1987 to 2002, most of the IS can be observed in Guangzhou and extended near the deviation of the PRD, where both the growth rate and speed of IS decreased, while after 2007, the IS expanded to cover the whole city, except the mountain and forest region. Furthermore, the growth rate and speed of the IS continued to increase, then decreased a little in 2017.

It can be observed that in Dongguan, IS experienced expansion based on multicore in the early period due to the policy of “Rural Economy.” The growth rate and speed of the IS experienced a “down-up-down” trend, while after 2002, the growth speed became faster and the urban construction was mainly around the core region, with the region closed to Shenzhen expanding rapidly. The urban expansion of Huizhou, Jiangmen, Zhaoqing, and Zhongshan presented a trend of the existing core regions driving the construction and expansion of new town. In Huizhou, the growth rate of the IS experienced a “down-up-down” trend, and the growth speed showed “W”-type development. The IS increasingly spread from the middle and southwest region to the western part of the city during 1987-2002. After 2002, northwest expansion can be observed, and the growth speed increased significantly. However, the growth rate and speed of IS in Jiangmen fluctuated with a high margin of growth speed and the IS clusters were dispersed during 1987-1997. Furthermore, the IS concentrated in the northeastern region and then expanded to the whole city. In Zhaoqing, the growth rate and speed increased firstly and then decreased, and the growth speed after 2002 was significantly faster. However, the increase in IS concentrated in the northwestern and southeastern regions of the city, while the growth rate of the IS in Zhongshan tends to decrease gradually but growth speed experienced an “up-down” trend.

4.2. The Variations of IS in Guangdong-Hong Kong-Macao Greater Bay Area Based on the Multispace Scale
4.2.1. The Scale of Political Institution

Based on the PRD integrated development zone, Hong Kong and Macao were included in the urban agglomeration of GBA. However, due to the different administrative nature, the IS change differs, as shown in Table 3. From 1987 to 2017, the IS of the PRD increased with an average annual growth speed of 345.91 km2/a with massive expansion. Moreover, the growth speed of IS from 2002 to 2017 is faster, with a growth speed of 616.56 km2/a of IS from 2007 to 2012.

The nine cities in the PRD have completed the transformation from the traditional agricultural economic regions to the important manufacturing centers. The IS of Hong Kong and Macao increased with an annual average speed of 3.72 km2/a, and they expanded mainly on the basis of the old city due to the limitation of administrative size; the growth speed of IS from 2002-2017 is slightly higher than that from 1987-2002. Moreover, The IS experienced rapid growth with a speed of 8.21 km2/a during 2012-2017. In the beginning, the urbanization processes are totally different between the PRD and special administration regions due to the different political systems. However, after implementation of “Closer Economic Partnership Arrangement (CEPA)” [80] policy in 2003, “Outline of the reform and development plan for the pearl river delta region (2008-2020)” policy in 2008 [81], and the framework for Guangdong-Hong Kong and Guangdong-Macao cooperation, the overall IS structure changed dramatically, and the variation tendencies of the growth rate of IS in the PRD and special administration regions were similar. However, the highest growth rate of increasing IS can be witnessed during 1987-1992.

Coordinated development of urban agglomeration has become the major impetus of global economic. Development and construction of GBA have become the national development strategy of China. We have found that there is a huge difference of urban evolution between the PRD and special administration regions. Constructing and developing GBA based on the differentiation of political institution and economic foundation have a significant impact on China’s economic development.

4.2.2. The Scale of Metropolitan Areas

In mainland China, around the GBA, there are three metropolitan areas: (1) Guang-Fo-Zhao: Guangzhou, Foshan, and Zhaoqing; (2) Shen-Guan-Hui: Shenzhen, Dongguan, and Huizhou; and (3) Zhu-Zhong-Jiang: Zhuhai, Zhongshan, and Jiangmen. To reveal the spatiotemporal characteristics of IS, the increasing area and the growth speed from 1987 to 2017 are calculated and presented in Figure 8. The IS of Guang-Fo-Zhao from 1987 to 2012 experienced rapid development; then the relative slow-down trend of IS growth can be seen after 2012; the IS area of Guang-Fo-Zhao increased 4215.13 km2. The IS area of Shen-Guan-Hui increased 2921.52 km2 over time. However, there are two peaks of growth speed from 1987 to 1992 and from 2002 to 2007, but the growth speed of IS decreased sharply after 1987-1992. Moreover, the growth speed experienced a slow-down trend followed by an increasing trend after 2002-2007. The IS of Zhu-Zhong-Jiang increased 3349.76 km2, the growth speed increased with an “up-down-up-down” trend, and the fastest growth speed can be observed in 2007-2012. However, the IS areas of these three metropolitan areas increased over time; the highest range of increment and fastest growth of IS areas occurred during 2007-2012 in Guang-Fo-Zhao with an increased speed of 298.16 km2/a. Moreover, Zhu-Zhong-Jiang experienced a period of rapid urbanization from 2007 to 2012 with an increased speed of 235.51 km2/a. Additionally, for Shen-Guan-Hui, the rapid urban growth can be observed during 1987-1992 and 2002-2007, with an increased speed of 138.58 km2/a and 151.64 km2/a.

To further explore the spatiotemporal change of IS for these three metropolitan areas, the SDE and GC of IS during 1987-2017 are calculated and are shown in Figure 9. The obvious directivity of the IS over time can be noticed. (1)In Guang-Fo-Zhao, from 1987 to 2002, the azimuth is distributed at 110°. However, from 2002 to 2007, the azimuth shifted to 114.04° with an increased speed of 0.85°/a, which indicates that the IS increased with a directivity trend towards northwest-southeast orientation. Moreover, the azimuth maintained around 112° after 2012. Furthermore, the flattering ratio first increased and then decreased, indicating that the IS area is distributed mainly along the long axis. There was a weak polarization phenomenon towards the short axis after 2012. The areas of ellipse deviation increased volatility over time; this indicated that the imperious surface continued to expand outward. It can also be observed that the IS GC moved from the border between Guangzhou and Foshan towards the border between Foshan and Zhaoqing in the past 30 years, which is mainly close to the geometric center. Thus, the GC moved 31.73 km northwest. This trend also indicates that the IS of Guang-Fo-Zhao increased mainly in the downtown of Guangzhou then extended towards the border between Guangzhou and Foshan in the 1990s. Finally, after the implementation of the Guang-Fo city plan in 2002 [82], the IS in Foshan expanded in the whole city rapidly, which created an integrated spatial structure of Guangzhou and Foshan. Meanwhile, the IS surface concentration in Zhaoqing started to undertake the industrial and labor transfer from Guangzhou and Foshan, thus indicating that, soon, Guangzhou, Foshan, and Zhaoqing will be integrated(2)In Shen-Guan-Hui from 1987 to 2002, the azimuth shifted and thus decreased from 73.16° to 70.74°. There is a significant counter-clockwise shift, where the azimuth decreased to 55.74° in 2007. This indicates that the orientation of IS distribution changed from east-west to northeast-southwest. However, the azimuth is maintained at 72° from 2012 to 2017, which shows that the spatial direction of IS expansion is towards east-west again, but this trend is gradually weakening. Moreover, the flattering ratio shows an M-type decreasing trend from 1.21 in 1987 to 1.15 in 2017, which reflects an uncertain principle expansion direction of the IS. Furthermore, the IS GC laid between Huizhou and Dongguan border and moved 3.16 km towards the northeast. In order to overcome the limitation of available space in Shenzhen, the Shen-Guan-Hui metropolitan area was established so that Dongguan and Huizhou can provide developed space to Shenzhen [83]. However, the cooperative mechanism was loose because of the distance between the downtowns in these three cities. Thus, the spatiotemporal pattern changing characteristic of IS was not obvious. Shen-Guan-Hui promoted the development of regional integration by reconstructing the spatial structure(3)In Zhu-Zhong-Jiang, the azimuth experienced an increasing-decreasing trend, and the azimuth shifted to 76.93° at an increasing speed of 0.64°/a from 1987 to 2007, which indicated that the IS expanded towards northeast-southwest orientation. However, this spatial structure started to decrease with 73.55° azimuth in 2017. Moreover, the flattering ratio presents a decreasing-increasing-decreasing trend from 1.52 to 1.82. Furthermore, the GC for Zhu-Zhong-Jiang is concentrated in Jiangmen, which surrounded the geometric center in U-type distribution, and the GC moved to the northeast with the movement of 13.21 km. There is no core city in Zhu-Zhong-Jiang and is only considered a geographical combination. The orientation of IS distribution moved from the northeast-southwest with the development of Macao and Zhuhai. Also, IS expansion is affected by the geomorphology and environment of the coastal region.

Generally, the spatial and temporal patterns of metropolitan areas in GBA are different. The IS distribution of Guang-Fo-Zhao and Shen-Guan-Hui has a high degree of aggregation, but IS distribution of Zhu-Zhong-Jiang is relatively dispersed. Zhu-Zhong-Jiang has natural and ecological advantage among these three metropolitan areas. The economic scale and industrial innovation strength of Guang-Fo-Zhao and Shen-Guan-Hui are in the domestic leading position. Multicore and group developing weaken the boundaries of cities in GBA, and the boundaries between downtowns and rural regions in GBA have a fusion trend. In the future, these three metropolitans should take industrial development and environmental protection into account, to solve the uneven resource allocation and urban development, and ultimately achieve the integrative development of GBA.

4.2.3. The Scale of the Development Poles

In order to improve the integral strength and global influence of GBA, three development poles (namely, Guangzhou-Foshan, Hong Kong-Shenzhen, and Macao-Zhuhai) were created for the GBA [61]. In the current study, the expansion area and growth speed of the development poles are calculated and presented in Figure 10. The IS of Guangzhou-Foshan increased by 2926.91 km2 during 1987-2017, with an average annual growth speed of 97.56 km2/a.

The SDE and GC of IS on three development poles during 1987-2017 are calculated and shown in Figure 11. (1)During 1987-2017, the azimuth decreased from 52.76° to 51.17°, which indicates that the IS of Guangzhou-Foshan is distributed towards the northeast-southwest. The flattering ratio increased during 1987-2017, which indicated that the IS expanded towards the long axis. The area of SDE increased overtime, revealing that the IS experienced rapid development. The GC of the IS was mainly situated towards the border of Guangzhou and Foshan, and the GC moved 5.29 km towards the southeast. However, near the downtown of Guangzhou and Foshan, the increasing IS is distributed. The spatial structure of Guangzhou-Foshan urban integration became stable after 2012; then, the IS started to extend towards the coastline(2)In Hong Kong-Macao, the azimuth shifted from 125.07° to 85.69° with a decreasing angle of 39°, which revealed that the orientation of the distribution of the IS during 1987-1992 shifted from northwest-southeast to east-west. However, the azimuth increased to 115.69° in 1997, which revealed that the IS shifted towards northwest-southeast. Moreover, the azimuth experienced an increasing trend indicating that the northwest-southeast spatial structure was enhanced during 2002-2017. The flattering ratio decreased from 1.27 to 1.10 overtime; in other words, the short axis length was increasing closer to the long axis length, and the area of SDE increased over time, which reflected the increasing diffusion degree of IS. The GC of Hong Kong-Macao is mainly distributed in Shenzhen towards the north of the geometric center of the development pole and moved towards the southeast with a distance of 3.30 km. In general, the increased IS was mainly distributed along the border of Hong Kong and Shenzhen; the IS of Shenzhen experienced a multicenter and leap forward development because of the flat terrain. However, Hong Kong, as a multi-island city, experienced an expansion by leaps and bounds without the obvious directivity, thus making it hard to understand the increase in the directivity(3)The azimuth of Macao-Zhuhai decreased from 104.94° to 95.26°, which exhibits the weak spatial structure of IS towards the northwest-southeast. The flattering ratio decreased from 1.84 to 1.56 during 1987-2007. Moreover, it increased to 1.68 in 2012, and then it decreased slightly in 2017, which shows the change in the SDE area with obvious directivity of the distribution of the IS. The GC of Macao-Zhuhai is located in Zhuhai, while the geometric center is located west of Macao-Zhuhai. However, the GC moved 4.42 km towards the southwest. Macao-Zhuhai and Hong Kong-Shenzhen have a similar trend of expansion, while in the beginning, the increasing IS towards the border of Macao and Zhuhai can be observed. Furthermore, the IS extended towards the whole region of Zhuhai

Currently, the development of GBA has formed a pattern that development poles promote the development of whole region. The cities in development poles of GBA are well integrated that prompted the regional cooperative construction and industrial economic development. Although integrative development of GBA is taking shape due to the positive influence of development poles, the connections of development poles and other cities in GBA are still weak because of their great disparity of development. Therefore, development poles should continue to lead industry and play important roles on the development of the whole region.

4.2.4. The Scale of the Core Cities

Hong Kong, Macao, Guangzhou, and Shenzhen are considered to be the “core engines” [61] of GBA. However, the process of change and the development pattern are different due to different scale and nature. Compared with that of Guangzhou and Shenzhen, the trend of IS expansion of Hong Kong and Macao was relatively stable and slow due to the rapid urbanization. To better analyze and understand the change of urban internal spatial structure of the “core engines,” the Patch Density (PD), Landscape Shape Index (LSI), and Shannon’s Diversity Index (SHDI) of the landscape pattern indexes are utilized in the current study. The results of the landscape pattern index of IS in the core cities of GBA from 1987-2017 are shown in Figure 12. (1)The PD and LSI of IS in Guangzhou decreased with the increase of IS density, which implies that the higher the density of the IS, the stronger is the aggregation. It also indicates a simple and more stable structure of the IS. The PD of very low-density IS decreased, and others did not change a lot. The LSI of very low-density IS decreased then increased while others increased. The SHDI of IS in Guangzhou experienced an M-type trend, which indicated that the IS in Guangzhou experienced the replacement process of the landscape. Finally, the natural surface had a high dominant position, but it changed to high-density IS(2)From 1987 to 2017, the PD and LSI of low-density IS in Hong Kong experienced a “decreasing-increasing” trend, but the very low-density IS had an opposite trend, and the trend of other densities is not steady, indicating that the degree of fragmentation and aggregation of different densities IS are not stable. Moreover, the PD of Hong Kong was lower than the PD of Guangzhou and Shenzhen, but LSI values were similar to those of Shenzhen, which indicates that the IS of Hong Kong had a lower degree of fragmentation. The SHDI of Hong Kong showed a similar changing trend with that of Guangzhou before 2002. However, the SHDI decreased again during 2007-2017, which might resulted in the replacement of lower density IS with higher IS(3)In Macao, the PD and LSI of lower density IS are higher than those of higher density IS, the PD of very low-density IS had a similar variation trend with Guangzhou, and the PD of other densities IS varied, indicating that the degree of fragmentation of IS declined gradually. The LSI of very low and low density IS increased first, then decreased and finally reached a stable state, which indicates that the degree of aggregation increased first then decreased; the LSI of other densities’ IS had a fluctuating trend with low values. Moreover, the value of PD and LSI is lower than those of Guangzhou and Shenzhen, which indicates that the IS had a low degree of fragmentation and a high degree of aggregation. Due to the “up-down-up” trend, the SHDI of Macao started a steady development; it shows that in 1987, the different densities’ IS had a uniform pattern. Finally, it can be observed that the higher IS had a high dominant position in 1992 with a homogenized pattern of IS(4)Shenzhen had a similar PD variation trend with Guangzhou, and PD of very low-density IS decreased from 1987 to 2017; others did not change a lot; it means that the degree of fragmentation of IS declined gradually. LSI of very low-density IS experienced the process of “decreasing-increasing”; others went through an “increasing-decreasing-increasing” process while the change range was smaller than that of Guangzhou, indicating that the degree of aggregation tends to be stable. The change characteristic of SHDI in Shenzhen was similar to that in Guangzhou. However, it tended to be stable after 2002, which indicated that the IS pattern is in a stable state

In summary, both Guangzhou and Shenzhen experienced a replacement process of the main landscape. Therefore, the dominant position of the spatial structure changed from lower IS to higher IS. The landscape of Hong Kong and Macao had a relatively low degree of fragmentation, and a high degree of aggregation compared with those of Guangzhou and Shenzhen and the spatial structure of the IS is more stable. In the future, Guangzhou, Shenzhen, Hong Kong, and Macao should pay more attention to rational urban development and ecological environment protection and also focus on optimizing the industrial structure to further enhance their effects for other regions in GBA. With the positive influence of core cities, GBA will become a real world-class bay area.

5. Conclusions

Information on IS area change in the context of rapid urbanization is key to understanding the process of urban agglomeration and deal with these serious challenges in terms of the environment, climate, and natural resources. However, previous studies have mostly focused on spatial structure changing in the city scale. In this study, we focused on GBA and investigated how the IS expanded during 1987, 1992, 1997, 2002, 2007, 2012, and 2017 at regional and city scales by using sensor-based Landsat images (i.e., TM, ETM+, and OLI). Furthermore, we discussed the spatiotemporal change characteristics of GBA growth at different spatial scales and indicated the evolution process of IS. Moreover, we attempted to understand the spatial structure evolution of urban agglomeration. We found that GBA has experienced a spatial structure evolution as follows: mega city, metropolitan area, and an urban agglomeration with rapid urbanization and industrialization. In addition, GBA has experienced a spatial reconstructing during its growth process, which is different compared with other urban agglomeration.

The results indicate that the IS of GBA experienced a rising trend during the past 30 years and exhibited a distinct characteristic along the PRD and the coastline. The IS area of GBA increased 749.68% during 1987 to 2017, the IS distributed along the Pearl River and coastline, and the ratio of the impervious surface decreased with the increasing distance from the axis. Furthermore, the IS of different cities in GBA changed. However, small cities developed by leapfrog development, while developed cities tended to expand on the edge of exiting IS. The findings of the current study reveal the evolution process in urban agglomeration. It not only provides the reference for GBA planning management but also supports a reference for a better understanding of other urban agglomerations. The current study only focused on multiscale analysis to observe the spatiotemporal patterns of the IS of urban agglomeration and preliminarily discussed the causes. However, the driving factors of sustainable development and ecological influence of urban agglomeration are more complicated with different political institutions, functions, and degrees of development.

These change progressions and orientations in IS patterns during the study period show the urbanization evolution from the pattern to the process in GBA, which thus provides useful information for urban sustainable planning and environmental protection. Further studies may be conducted on revealing them quantitatively.

Data Availability

The topographic map of the Institute of Surveying and Mapping, Department of Natural Resources of Guangdong Province and the Digital Orthophoto Map (DOM) data of Guangzhou Jiantong Surveying, Mapping, and Geoinfo Co., Ltd., used to support the findings of this study were supplied by Xiaoding Liu and Mingqiang Chen under license and so cannot be made freely available. Requests for access to these data should be made to Fan Liu ([email protected]).

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The study was supported by the National Natural Science Foundation of China (41871292), the Science and Technology Program of Guangdong Province, China (2018B020207002), and the Science and Technology Program of Guangzhou, China (201803030034 and 201802030008).