Impervious Surface Expansion: A Key Indicator for Environment and Urban Agglomeration—A Case Study of Guangdong-Hong Kong-Macao Greater Bay Area by Using Landsat Data

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.


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 [2][3][4][5]. 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" [7][8][9]. 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 [14][15][16][17][18] 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 [23][24][25]. 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 [32][33][34].
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., multitem-porality, 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 [47][48][49][50][51]. 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 crossadministrative 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).

Journal of Sensors
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.

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 "3 + 3 + 3" 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 km 2 , 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.   [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.

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 1987GBA in , 1992GBA in , 1997GBA in , 2002GBA in , 2007GBA in , 2012GBA in , 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. 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.    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: Overall accuracy = TP + TN TP + TN + FP + FN , Sensitivity = TP TP + FN , 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.

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 "σ x " and "σ y " in directions "x" and "y" are calculated as follows: where ðx i , y i Þ is the spatial location of the object, w i is the corresponding weight, ðx i , y i Þ is the coordinates deviation from the spatial location, and ( X w , Y w ) is the gravity center.

Journal of Sensors
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.  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 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].  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 21 st 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.

Results and Discussion
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.  Representing the complexity of spatial structure, the higher PD is, the more fragmental the class is [76].
In the measurement of patch aggregation, the larger LSI is, the more discrete the patches are [75].
Landscape level SHDI 0 ≤ SHDI Representing the richness and complexity of the land surface, the higher the value, the higher the complexity of the whole landscape [77]. 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.  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   Journal of Sensors 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.  [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.

The Variations of IS in
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.     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 out-ward. 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°a zimuth 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   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.

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 km 2 during 1987-2017, with an average annual growth speed of 97.56 km 2 /a.
The SDE and GC of IS on three development poles during 1987-2017 are calculated and shown in Figure 11.         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.

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  (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 "updown-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.

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.