Co-Evolution of Emerging Multi-Cities: Rates, Patterns and Driving Policies Revealed by Continuous Change Detection and Classiﬁcation of Landsat Data

: The co-evolution of multi-cities has emerged as the primary form of urbanization in China in recent years. However, the processes, patterns, and coordination are not well characterized and understood, which hinders the understanding of the driving forces, consequences, and management of polycentric urban development. We used the Continuous Change Detection and Classiﬁcation (CCDC) algorithm to integrate all available Landsat 5, 7, and 8 images and map annual land use and land cover (LULC) from 2001 to 2017 in the Chang–Zhu–Tan urban agglomeration (CZTUA), a typical urban agglomeration in China. Results showed that the impervious surface in the study area expanded by 371 km 2 with an annual growth rate of 2.25%, primarily at the cost of cropland (169 km 2 ) and forest (206 km 2 ) during the study period. Urban growth has evolved from inﬁlling being the dominant type in the earlier period to mainly edge-expansion and leapfrogging in the core cities, and from no dominant type to mainly leapfrogging in the satellite cities. The unfolding of the “cool center and hot edge” urban growth pattern in CZTUA, characterized by higher expansion rates in the peripheral than in the core cities, may signify a new form of the co-evolution of multi-cities in the process of urbanization. Detailed urban management and planning policies in CZTUA were analyzed. The co-evolution of multi-cities principles need to be studied in more extensive regions, which could help policymakers to promote sustainable and livable development in the future.

. The location of the study area: (a) the study area in China; (b) the study area in the Hunan Province; (c) the distribution of the nine regions in Chang-Zhu-Tan.

Data and Pre-Processing
In this study, all available Level 1 Terrain (L1T) Landsat 5, 7, and 8 images for the World Reference System 2 (WRS2) path/row 123/040 and 123/041, with the percentage of cloud cover below 80%, were obtained ( Figure A1). There was a total of 637 images meeting these criteria from 2001 to 2017. The image preprocessing was completed using ESPA (https://espa.cr.usgs.gov/), which provides an order mode to batch processing Landsat images [62]. All images were atmospherically corrected to the surface reflectance with established algorithms and presented a high consistency between the different sensors [63][64][65]. Specifically, Landsat 5 and 7 surface reflectances were processed by the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm (version 3.4.0) [66] and the Landsat 8 surface reflectances were processed by the Land Surface Reflectance Code (LaSRC) (version 1.4.1) [67]; auxiliary input data, such as water vapor, ozone, atmospheric pressure, and aerosol optical thickness, were retrieved from Moderate Resolution Imaging Spectroradiometer (MODIS); and the digital elevation derived from the Earth Topography Five Minute Grid [68]. For each Landsat image, six surface reflectance bands (Blue, Green, Red, NIR, SWIR1, and SWIR2), one thermal band (Brightness Temperature) and one quality assessment band were employed. To detect clouds, cloud shadows, and snow, we used the Function of Mask algorithm (FMask) combined with potential cloud pixels and the cloud probability mask to derive the potential cloud layer through discarding most unusable observations to reduce noise in the analysis [69,70]. We downloaded the CCDC version 12.30 from the Global Environmental Remote Sensing (GERS) Laboratory, and the data preparations were finished by the CCDC workflow and supported by the software ArcGIS 10.2 and ENVI 5.1 [54]. The overview flowchart is shown in Figure 2. City and county boundary datasets were obtained from the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC) [22,71]. Auxiliary data were obtained from FROM-

Data and Pre-Processing
In this study, all available Level 1 Terrain (L1T) Landsat 5, 7, and 8 images for the World Reference System 2 (WRS2) path/row 123/040 and 123/041, with the percentage of cloud cover below 80%, were obtained ( Figure A1). There was a total of 637 images meeting these criteria from 2001 to 2017. The image preprocessing was completed using ESPA (https://espa.cr.usgs.gov/), which provides an order mode to batch processing Landsat images [62]. All images were atmospherically corrected to the surface reflectance with established algorithms and presented a high consistency between the different sensors [63][64][65]. Specifically, Landsat 5 and 7 surface reflectances were processed by the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm (version 3.4.0) [66] and the Landsat 8 surface reflectances were processed by the Land Surface Reflectance Code (LaSRC) (version 1.4.1) [67]; auxiliary input data, such as water vapor, ozone, atmospheric pressure, and aerosol optical thickness, were retrieved from Moderate Resolution Imaging Spectroradiometer (MODIS); and the digital elevation derived from the Earth Topography Five Minute Grid [68]. For each Landsat image, six surface reflectance bands (Blue, Green, Red, NIR, SWIR1, and SWIR2), one thermal band (Brightness Temperature) and one quality assessment band were employed. To detect clouds, cloud shadows, and snow, we used the Function of Mask algorithm (FMask) combined with potential cloud pixels and the cloud probability mask to derive the potential cloud layer through discarding most unusable observations to reduce noise in the analysis [69,70]. We downloaded the CCDC version 12.30 from the Global Environmental Remote Sensing (GERS) Laboratory, and the data preparations were finished by the CCDC workflow and supported by the software ArcGIS 10.2 and ENVI 5.1 [54]. The overview flowchart is shown in Figure 2. City and county boundary datasets were obtained from the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences GLC and included the ground truth samples for the classification and reference data for the temporalspatial accuracy assessment [72].

CCDC Algorithm and Annual LULC maps
Annual land cover changes were detected and mapped using the CCDC algorithm [37]. First, all available Landsat observations for each pixel were used to develop time series models with the seasonality, trend, and break components to detect the intra-annual, inter-annual, and abrupt changes by the following Equation (1) [37,73].
( , ) = 0, + 1, cos ( 2 ) + 1, sin ( 2 ) + 1, where F is the fitted model for band i and time x (Julian date), 0, is the coefficient for the mean of band i, 1, and 1, are the coefficients representing intra-annual change, 1, is the coefficient representing the inter-annual change, or trend, and T is the number of days per year (T = 365.25). If the land cover pattern does not change, the time series model will be continued. Once, the breaks were found by fitting a linear model to a stable history period of up to two dates. New observations were added and their residuals compared to the RMSE (Root Mean Square Error) of the history period. Whenever the error in the new observation was too high, the algorithm would place a break point and regard it as a potential change. If breaks accumulated enough in a point, a change was affirmed and a new stable model would be initialized. Generally, six breaks per pixel are used to provide a high confidence when the land-use pattern changed. The coefficients in Equation (1) were recorded as the input for the second part of the CCDC algorithm: land cover classification. LULC classification in the CCDC was accomplished using a Random Forest Classifier (RFC) [74], because of its robustness in mapping large-area land cover, especially effective when a large number of features are included. The RFC model parameters were the model coefficients and RMSE. To balance computation time and classification accuracy, 500 bagging iterations were implemented in the RFC procedure, and the number of variables to split each node was limited to the square root of the total number of explanatory variables. The output of the RFC for each year was a six-category land cover map to satisfy our demand. According to the FROM-GLC Level 1 Type land cover system [72] and the actual situation in CZTUA, the LULC in the study was classified as Cropland (FROM-GLC classes of Crop), Forest (FROM-GLC classes of Forest), Grassland (FROM-GLC classes of Grass,

CCDC Algorithm and Annual LULC maps
Annual land cover changes were detected and mapped using the CCDC algorithm [37]. First, all available Landsat observations for each pixel were used to develop time series models with the seasonality, trend, and break components to detect the intra-annual, inter-annual, and abrupt changes by the following Equation (1) [37,73].
where F is the fitted model for band i and time x (Julian date), a 0, i is the coefficient for the mean of band i, a 1, i and b 1, i are the coefficients representing intra-annual change, c 1, i x is the coefficient representing the inter-annual change, or trend, and T is the number of days per year (T = 365.25). If the land cover pattern does not change, the time series model will be continued. Once, the breaks were found by fitting a linear model to a stable history period of up to two dates. New observations were added and their residuals compared to the RMSE (Root Mean Square Error) of the history period. Whenever the error in the new observation was too high, the algorithm would place a break point and regard it as a potential change. If breaks accumulated enough in a point, a change was affirmed and a new stable model would be initialized. Generally, six breaks per pixel are used to provide a high confidence when the land-use pattern changed. The coefficients in Equation (1) were recorded as the input for the second part of the CCDC algorithm: land cover classification. LULC classification in the CCDC was accomplished using a Random Forest Classifier (RFC) [74], because of its robustness in mapping large-area land cover, especially effective when a large number of features are included. The RFC model parameters were the model coefficients and RMSE. To balance computation time and classification accuracy, 500 bagging iterations were implemented in the RFC procedure, and the number of variables to split each node was limited to the square root of the total number of explanatory variables. The output of the RFC for each year was a six-category land cover map to satisfy our demand. According to the FROM-GLC Level 1 Type land cover system [72] and the actual situation in CZTUA, the LULC in the study was classified as Cropland (FROM-GLC classes of Crop), Forest (FROM-GLC classes of Forest), Grassland (FROM-GLC classes of Grass, Shrub), Wetland (FROM-GLC classes of Wetland), Water (FROM-GLC classes of Water), and Impervious Surface (FROM-GLC classes of Impervious, Bareland). We set the dynamics from the impervious surface as the urban expansion indicator. The ground truth samples were collected via visual interpretation, aided by very high resolution (VHR) images from Google Earth TM and prior research [22,75]. We randomly selected 26,670 samples within the study area among the six land cover categories; the samples proportion for each LULC category was controlled by the percentage of each class area in FROM-GLC10 (2017) [54,76]. Dividing the samples into two groups, 90% of the samples were used to train the RFC, and the remaining 10% were used for the accuracy assessment. The user's and producer's accuracy were calculated for each land cover type along with the overall accuracy [77].

Annual Urban Expansion Rate
In order to study the characteristics of urban expansion quantitatively (e.g., speed, extent, intensity, and direction), the annual increment (AI) (km 2 year −1 ) and the annual urban growth rate (AGR) (%) were calculated and used to quantify and compare the urban expansion speed of CZTUA during the study period. AI could show the newly developed impervious surface area per year and AGR could be used to make a better comparison among the cities, which could eliminate the influence of different primary city sizes [10]. The two indices were defined as follows: where A start and A end are the urban area at the start and end year, respectively, while d is the time span of the study.

Urban Expansion Types
Urban growth patterns were generally classified as infilling, edge-expansion, and leapfrogging by Forman (1995) [78]. Infilling is defined as a new urban patch formed via filling in the gaps within existing urban patches. Edge-expansion is a type of urban growth in which a newly developed urban patch extends outward along the edge of existing urban patches. Leapfrogging growth represents new urban patches developed independently and without any existing urban patches. The landscape expansion index (LEI) [79] was calculated to define these three urban growth types using the following equation: where A 0 is the area of the intersection between the buffer zone of the new patch and the existing patches (occupied category) and A v is the area of the intersection between the buffer and the vacant category. Infilling growth patch is defined by an LEI larger than 50 and an edge-expansion growth patch is defined by an LEI smaller than 50 but not equal to zero. The patches with a zero LEI were defined as leapfrogging growth.

Landscape Metrics
In order to reveal the evolution of the landscape pattern in the process of urbanization, six prominent landscape metrics based on spatial and complexity criteria were calculated to illustrate the characteristics of the CZTUA, which could avoid the influences of the different scales and extents of the regional boundary. These metrics include the Percentage of Landscape (PLAND), Largest Patch Index (LPI), Landscape Shape Index (LSI), Patch Density (PD), Patch Cohesion Index (COHESION) and Interspersion and Juxtaposition Index (IJI). These metrics could reflect the impervious surface patches' spatial distribution and connectivity, which are relevant to urban expansion and able to show the results effectively [80]. PLAND describes the proportion of impervious surface in the total region area. LPI represents the percentage occupied by the largest patch, reflecting the city center. LSI and PD demonstrate the edge and density in the landscape, respectively. With impervious surface increase and the landscape becoming fragmented, the LSI and PD will increase, while the LSI and PD may decrease when the impervious surface patches merged [23]. COHESION explains the physical connectivity between the different impervious surface patches. IJI indicates the complexity of the patches around the impervious surface patches [81]. All of the landscape metrics were calculated for the impervious surface [20,23,27]. These metrics were computed for each of the selected cities at the landscape class with the help of FRAGSTATS 4.2. A detailed description for each metric is listed in Table 1. The formulas for the landscape metrics are listed in Equations (A1)-(A6). Table 1. Landscape metrics used in this study.

Landscape Metric Abbreviation Description
Percentage of Landscape PLAND The percentage the landscape of the corresponding patch type.
Patch Density PD The number of patches of per 100 ha.
Largest Patch Index LPI Proportion of total area occupied by the largest patch of a land cover type.

Landscape Shape Index LSI
A modified perimeter-area ratio of the form that measures the shape complexity of the whole landscape or a specific cover type.
Interspersion and Juxtaposition Index IJI The overall distribution and juxtaposition of various patch types Patch Cohesion Index COHESION Measuring the physical connectedness of the corresponding patch type.

Accuracy of the LULC Classification
For the entire study area, we selected five years of the LULC maps (i.e., 2001, 2005, 2010, 2015, and 2017) for our spatiotemporal accuracy assessment. Validation samples were the remaining 10% of the samples. The classification results achieved an overall accuracy of 90.44%-92.31% across these five years ( Table 2  The classification results of the cropland, forest, and water classes were more accurate, while those of grassland and wetland were relatively poor. The difference in classification accuracy of different land covers may be due to the fact that the number of training samples for the grasslands and wetlands were less than those of the other land-cover types, typical for rare land cover types [82,83]. Besides, seasonal variations may also lead to the related low classification accuracy. For example, the probability of a misclassification between a cropland and forest was still high, and the grasslands and wetlands were confused easily. Nevertheless, the small differences in overall accuracy among the different years suggest that the CCDC model performed well in mapping the time series of LULC. So, the results could meet the requirement for monitoring and analyzing the urban expansion process. Figure 3 shows the Landsat images and resultant land cover maps in 2001 and 2017, respectively, in a typical area.

Temporal LULC Change from 2001 to 2017
Based on the yearly classification maps, the land cover in CZTUA has significantly and rapidly changed in the past 17 years (Figure 4). In general, the impervious surface area increased while the cropland and forest decreased continuously, and the water area, wetland, and grassland remained relatively stable. The land cover was composed of cropland of 5297 km 2 , forest of 4355 km 2 , grassland of 12 km 2 , wetland of 152 km 2 , water of 344 km 2 , and impervious surface of 808 km 2 in CZTUA in the year of 2001 (Table A1). Overall, an increase of 371 km 2 in impervious surface, and a decrease of 169 km 2 and 206 km 2 in cropland and forest, respectively, were observed from 2001 to 2017 in the study area ( Figure 4). Most of the impervious surface was gained from both forest and cropland and the transformation from wetland, water, and grassland was not obvious (Table A2).
Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 26 area ( Figure 4). Most of the impervious surface was gained from both forest and cropland and the transformation from wetland, water, and grassland was not obvious (Table A2). The impervious surface area increased from 808 km 2 to 1179 km 2 over the past 17 years with the net increase rate of 45.90% (Table A1). There was large spatial variation in the area of impervious surface increase among subregions: from 211 to 246 km 2   The impervious surface area increased from 808 km 2 to 1179 km 2 over the past 17 years with the net increase rate of 45.90% (Table A1). There was large spatial variation in the area of impervious surface increase among subregions: from 211 to 246 km 2 (16.59%) for Changsha City; 145 to 227 km 2 (56.55%) for Changsha County; 126 to 217 km 2 (72.22%) for Wangcheng District; 82 to 106 km 2 (29.27%) for Zhuzhou City; 31 to 64 km 2 (106.45%) for Zhuzhou County; 63 to 94 km 2 (49.21%) for Liling City; 83 to 103 km 2 (24.10%) for Xiangtan City; 57 to 110 km 2 (92.98%) for Xiangtan County; and 9 to 13 km 2 (44.44%) for Shaoshan City (Table A4).  (Table A3).

Spatiotemporal Dynamics of Urban Expansion
CZTUA demonstrated great spatial variability in urban development from 2001 to 2017 ( Figure  6). Changsha City, Xiangtan City, and Zhuzhou City are the belt-shaped cities that were initially developed along the river. In recent years, urbanization in the southern part of Changsha City, the northeast of Xiangtan City and the northwest of Zhuzhou City has been active, showing a trend of triangular convergence (Figure 6 c). For Changsha County and Wangcheng District, the primary core has been constantly expanding outward, especially in the eastern part of Changsha County and western part of Wangcheng District (Figure 6 a, b). Meanwhile, the urban area of Zhuzhou County and Xiangtan County expanded rapidly after the year 2010, resulting in a doubling growth. Shaoshan City and Liling City were developmentally restrained.

Spatiotemporal Dynamics of Urban Expansion
CZTUA demonstrated great spatial variability in urban development from 2001 to 2017 ( Figure 6). Changsha City, Xiangtan City, and Zhuzhou City are the belt-shaped cities that were initially developed along the river. In recent years, urbanization in the southern part of Changsha City, the northeast of Xiangtan City and the northwest of Zhuzhou City has been active, showing a trend of triangular convergence (Figure 6c). For Changsha County and Wangcheng District, the primary core has been constantly expanding outward, especially in the eastern part of Changsha County and western part of Wangcheng District (Figure 6a,b). Meanwhile, the urban area of Zhuzhou County and Xiangtan County expanded rapidly after the year 2010, resulting in a doubling growth. Shaoshan City and Liling City were developmentally restrained. northeast of Xiangtan City and the northwest of Zhuzhou City has been active, showing a trend of triangular convergence (Figure 6 c). For Changsha County and Wangcheng District, the primary core has been constantly expanding outward, especially in the eastern part of Changsha County and western part of Wangcheng District (Figure 6 a, b). Meanwhile, the urban area of Zhuzhou County and Xiangtan County expanded rapidly after the year 2010, resulting in a doubling growth. Shaoshan City and Liling City were developmentally restrained.  Figure 7 shows the dynamics for the proportion of three urban growth types (i.e., leapfrogging, edge-expansion, and infilling) of the newly developed impervious surface patches in CZTUA during the study period. It could be seen that there was large spatial and temporal variability during the urban expansion processes across regions. Over the entire study period, infilling was the dominant growth type with a steady trend to decrease, possibly as the gaps became filled. On the other hand, leapfrogging had a marginal occurrence at the beginning of the time series, to later become the dominant urban growth type in 2016. Edge-expansion remained more stable throughout the time series. The proportion of infilling reduced steadily from about 80% in 2001 to 40% in 2016 in Changsha City, while the proportion of edge-expansion increased smoothly from about 20% to 40% and the share of leapfrogging increased from about 0% to 10% (Figure 7). In 2016, leapfrogging increased sharply in Changsha County and Wangcheng District. Most of the newly developed urban land in Changsha County concentrated around the existing urban area and is distributed with the newly developed highways and airport in the region (Figure 3b). Figure 8 shows the annual changes in landscape metrics in CZTUA. The PLAND and LPI increased slowly for all nine regions during whole period. Changes in LPI suggested the proportion of the largest urban patch (the existed impervious surface) in each region. Specifically, Changsha City presented the highest proportion in PLAND and LPI with a steady increased rate among the nine regions, while the rest of the regions showed a similar pattern in PLAND and LPI though the values were lower than those for Changsha City. The LPI in Xiangtan City, Shaoshan City, and Wangcheng District has greatly increased from 2007 to 2010.

Landscape Changes during Urban Expansion
series. The proportion of infilling reduced steadily from about 80% in 2001 to 40% in 2016 in Changsha City, while the proportion of edge-expansion increased smoothly from about 20% to 40% and the share of leapfrogging increased from about 0% to 10% (Figure 7). In 2016, leapfrogging increased sharply in Changsha County and Wangcheng District. Most of the newly developed urban land in Changsha County concentrated around the existing urban area and is distributed with the newly developed highways and airport in the region (Figure 3 b). The landscape metrics PD and LSI described the complexity and fragmentation for the impervious surface. The PDs in Wangcheng District, Changsha County, and Zhuzhou City were relatively higher than those in other regions, with a slow decline in the early stage and rapid growth in the latest period. The remaining regions had a relatively lower PD and could form a different group. Changsha City was different from all other regions by presenting a steady decrease in PD over the past 17 years. Contrarily, Liling City, Xiangtan County, and Zhuzhou County increased in the latest period and had no significant changes in the earlier stage. In the total area of CZTUA, PD changed a little during the observation and kept a relatively low value. For the LSI, CZTUA had the highest LSI with the U-shape function changed. In this study, the LSI was lower than 150 in the past 17 years. According to the value of the LSI, the cities in CZTUA could be roughly divided into three groups, the first group of which were greater than 90, including Changsha County, Shaoshan City, Liling City, and Xiangtan County; for the second group, which includes Zhuzhou County, the LSI was between 60 to 90; the last group, with an LSI lower than 60, includes Zhuzhou City, Changsha City, Xiangtan City, and Shaoshan City.
Most of regions in CZTUA were under a steady and similar pattern in LSI. The IJI and COHESION metrics represent aggregation and connectivity. IJI in Changsha City was the highest throughout the period, rising from 50 in 2001 to just under 70 until 2013, when there was a levelling off for three years. The remaining regions presented a similar pattern, which decreased first and then increased or kept the value without change for IJI. Wangcheng District and Changsha County had a reverse U-shape during this observation. For COHESION, Changsha City, Xiangtan City, Zhuzhou City, Wangcheng District, and Changsha County had a relatively higher value, which was close to 100, and presented a steady value or slight increase from 2001 to 2017. Xiangtan County, Shaoshan City, and Liling City increased steady before 2009 and stabilized in the latest period, while Xiangtan County, with a little bit of hysteresis, persisted by increasing up to 2013. Zhuzhou County presents sustainable growth in COHESION.  Figure 8 shows the annual changes in landscape metrics in CZTUA. The PLAND and LPI increased slowly for all nine regions during whole period. Changes in LPI suggested the proportion of the largest urban patch (the existed impervious surface) in each region. Specifically, Changsha City presented the highest proportion in PLAND and LPI with a steady increased rate among the nine regions, while the rest of the regions showed a similar pattern in PLAND and LPI though the values were lower than those for Changsha City. The LPI in Xiangtan City, Shaoshan City, and Wangcheng District has greatly increased from 2007 to 2010. The landscape metrics PD and LSI described the complexity and fragmentation for the impervious surface. The PDs in Wangcheng District, Changsha County, and Zhuzhou City were relatively higher than those in other regions, with a slow decline in the early stage and rapid growth in the latest period. The remaining regions had a relatively lower PD and could form a different group. Changsha City was different from all other regions by presenting a steady decrease in PD over the past 17 years. Contrarily, Liling City, Xiangtan County, and Zhuzhou County increased in the latest

Co-Evolution of Multi-Cities in CZTUA
For the nine regions in CZTUA, Changsha County, Wangcheng District, and Xiangtan County had the greatest magnitude of urban expansion from 2001 to 2017, with the net increase area of 82 km 2 , Remote Sens. 2020, 12, 2905 14 of 24 91 km 2 , and 53 km 2 , respectively. Xiangtan county, Zhuzhou county, and Wangcheng District had the highest net increase rates, with 92.79%, 105.10%, and 71.91%, respectively (Table A4). Besides, Changsha City had the greatest primary impervious surface area, but its net increase rate was the smallest among the study area (Table A4).
Changsha City is the provincial city of Hunan province, which developed earlier with prior policy support. Cheap labor promoted manufacturing and heavy industry, and thousands of natural resources were developed in that early stage. The large existence of backward primary industry presented a great challenge for sustainable development pattern of the city [43]. Consequently, a strategy has been put in place to promote industrial upgrading, translocation, and co-evolution with neighboring cities [84]. Changsha County and Wangcheng District became the main receiving area for some of the city's industries, especially in 2003, 2007, and 2013. In these periods, Changsha City had slowed down the pace of urban expansion while Changsha County and Wangcheng District accelerated the pace of urban expansion and achieved the greatest AIs in CZTUA ( Figure 5). Both PD and LSI in Changsha City, Wangcheng District, and Changsha County decreased from 2001 to 2013, indicating that the impervious surface patches gradually merged (Figure 8). Besides, the COHESION of Zhuzhou county, Xiangtan City, and Xiangtan County increased beyond 3%, especially after 2009, indicating the connectivity in these regions achieved a high level in a short time (Figure 8). The AGR in Xiangtan County and Zhuzhou county had presented the highest values, reflecting that these regions were also promoted ( Figure 5). The co-evolution pattern of the multi-cities could be summarized as follows: the core city developed first due to policy advantages, and the surrounding cities developed primarily benefited from the industrial upgrading of the core city. As a result, all cities were developed and presented co-evolution pattern where the cities in a specific region play a coordinated urban function. Many social, economic, and industrial studies show a similar pattern [12,15,16,85].
The co-evolution pattern of multi-cities in CZTUA may be an example of city development in other regions as some international megalopolises experienced a similar urbanization process. For instance, the Northeast megalopolis, the most populous megalopolis in US, has spent a very long time in integrating the major cities such as Boston, New York City, and Washington, D.C., along with their metropolitan areas and suburbs [86]. The urban industrial upgrading and transformation from major cities in the later 20th century had promoted the development of the surrounding cities. The Atlantic Ocean coastal cities connected with each other from northeast to southwest and differentiated into diverse specific functions, which could be regarded as the co-evolution model [84,87]. China and other developing countries are following this urbanization pattern and promoting urban area integration, sustainability, and livability in relatively short term [42,88].

Spatiotemporal Dynamics of Urban Expansion and the Possible Drivers
CZTUA presents a triangle pattern in geography, which made the urban agglomeration develop a polycentric pattern. Besides, the Xiang River winds its way through the core cities and provides a natural urban developing pattern [61]. On one hand, urbanization brings vitality to the city, promoting residents' lives more efficiently and conveniently; but on the other hand, the negative consequences, especially in the old town, also emerged with urban expansion. Therefore, managing and promoting the urban areas to be livable and sustainable is a crucial and thorny issue for the management of the national and local governments. All these characteristics lead to distinct specific dynamics in the landscape metrics and urban expansion types during the urbanization process in CZTUA. In the earlier period, the three core cities were relatively independent, physically along the Xiang River, and the suburbs were not associated. As the transportation systems, such as highways, the G-series high-speed train, and urban rail, were built up, all regions in CZTUA were gradually connected, resulting in an increased IJI and COHESION.
The spatiotemporal difference between urban and suburban development was an important driving force of landscape change in CZTUA [89]. For urban areas, such as Changsha City, Zhuzhou City, and Xiangtan City, the dominant urban growth type had experienced an obvious transformation from infilling to edge-expansion. That would be explained by urban development theory [90,91], building main roads inside a city firstly and improving the connectivity between the urban and suburban areas, which also promote the inner-city grid gradually (Figure 6 a, b). Most newly developed impervious surfaces were distributed between the roads and filled-up blanks, resulting from the dominant urban growth and infilling (Figure 7). Subsequently, the inside city was gradually filled up, and with the continuous increase in land prices, the urban expansion type turned to edge-expansion, as shown by previous similar studies [10,23,27]. For the suburban areas, urbanization was relatively slow and natural, and there was no significant dominant urban growth type presented in the earlier period. With the three core cities built well, the expansion type was converted to edge-expansion because the capital and workforce transferred from urban to suburban [20]. In the latest period, a series high and new technology industrial development zones were designated in urban-rural fringe areas far away from the existing impervious surface (Figure 6c,d). These zones played a crucial role in both economic development and urban construction [92]. As a consequence, the urban expansion type shifted from edge-expansion to leapfrogging [61].

Urban Expansion under the Guidance of Policy
The nine regions in CZTUA had a steady urban increase over the past 17 years.  Figure 6, Table A4). These rates were highly dependent on the administrative level, which can be roughly divided into three types of growth: the stable-growing type, including Changsha City, Zhuzhou City, and Xiangtan City; the fast-growing type, including Wangcheng District and Changsha County; and the slow-growing type, including Shaoshan City and Liling City. These characteristics of the CZTUA had a very important and close overlap with national development strategies and economic policies, as well as the local government measures and implementations [93].
Earlier national policies for urbanization in China mainly focused on coastal megacities and land-locked cities did not catch up with the pace of development, therefore some cities in CZTUA were relatively backward. Until the year 2002, the Rise of Central China policy was adopted by the central government, mainly aimed at improving the central China's capacity for the coordinated development [2,94]. The local government of Hunan province, promoting the integration of CZTUA under a scientific and systematic urban planning, shaped up the first urban agglomeration plan in 2005. This plan considered the differences in physical geography and economic background of the three core cities, spelled out different urban designations and orientations, eventually made a significant impact on the evolution of the cities. Specifically, Changsha City is a comprehensive metropolis for the development of culture, education, finance, science, and technology, and played a leading role; Zhuzhou City focused on the development of producer services on the basis of transportation and industries; Xiangtan City was based on the heavy industry and has gradually promoted the incubation of high-tech industries, which is closer to the producer services. Three core cities along the Xiang River took the leadership in forming a megacity, while other regions were less driven by policies during this period. The two policies implemented in 2000 and 2005 promoted the growth of urban expansion from 2000 to 2010, of which Changsha City, Wangcheng District, and Changsha County contributed the largest proportion of growth, followed by Xiangtan City and Zhuzhou City, and the rest still grew slowly [95].
Later on, the resource-saving and environment friendly urban development strategy was proposed by both the central government and the local government in 2007. According to the different characteristics of the cities, different urban development goals were defined so as to achieve the suitable and sustainable development goals. It was suggested that the original policies should be continued in the three core cities and in the suburban areas, such as Shaoshan City and Liling City; also, a core conservation area connecting the three core cities has been designated as the green heart and any urban land development would be banned there (Figure 6c). Some prior urban developing regions, such as Wangcheng District, Changsha County, and Xiangtan County, have been designated primarily for economic development. In 2010, the National Economic and Technical Development Zones were approved and setup in Jiuhua, Xiangtan County, resulted in large-scale expansion of the city ( Figure 5) [96]. Since 2012, the ecological red lines policy has been implemented in CZTUA, which curbed illegal and activities (e.g., mining in the core conservation regions) and even converted some urban land (e.g., brown fields) to vegetated surfaces, which leading to negative numbers during the years from 2013 to 2015 in Changsha City ( Figure 5). Managing urban development in accordance with location-specific features has become the main ideology of the current policy.
The significant difference between the core cities and the surrounding cities may be caused by the guidance of the national and local policies. The observations showed that the core cities have grown steadily over the past 17 years, indicating the urban expansion rates have slowed down as these regions approach their peak in urbanization process. In the prior developing satellite cities, social capital from the core cities has been attracted to accelerate their urbanization process. We concluded this pattern as "cool center and hot edge", describing the pace of the urban expansion in the core cities that has slowed but having been accelerated in the satellite cities. This pattern may indicate a transformation in co-evolution of multi-cities with powerful policy guidance in China.

Conclusions
Urbanization has been one of the most significant human impacts on the environment, and it is believed to be one of the biggest challenges for the national and local government in China. A coordinated development pattern in polycentric urban agglomeration is regarded as a solution to avoid the adverse impacts from urbanization. However, there are few studies that could reveal the spatiotemporal details of the urbanization, resulting from inferior data and inappropriate data processing tools. This study provided a continuous and comprehensive understanding of the spatiotemporal landscape dynamics in a typical urban agglomeration of China, which was timely and necessary to support the regional urban master planning.
This study used the CCDC algorithm to integrate all available Landsat images and obtain annual LULC maps. Since this algorithm was pixel-based, it required substantial computing power but could effectively capture the annual landscape dynamics in the urban area. Tracking the LULC changes could provide a unique perspective for decisionmakers and promote more sustainable, livable urban development. The results showed that CZTUA has experienced a rapid change over the past 17 years. A total 371 km 2 (with annual growth rate 2.25%) impervious surface expansion was observed with cropland and forest decreasing by 169 km 2 and 206 km 2 in the CZTUA, respectively; most of the impervious surface was converted from cropland and forest. Annual growth rates slowed down gradually in the core cities and accelerated in the surrounding areas; we concluded this pattern as "cool center and hot edge", which may indicate an important urbanization transformation pattern in China. Urban expansion types also had a significant transformation. The dominant urban expansion type changed from infilling to leapfrogging in the three core cities. While not significant in the earlier period, edge-expansion and leapfrogging became the main urban expansion types in the latest period in the suburban areas. The co-evolution pattern funded in this study indicated that multi-cities may experience an inevitable stage in the process of urbanization and more further studies need to be conducted to explore and summarize more extensively in different countries with different policy backgrounds and in the different stages of urban development. Formulas of the Landscape metrics in this study.
where P i is the proportion of the landscape occupied by patch type i, a ij is the area (m 2 ) of patch ij, A is the total area (m 2 ) (A1)(A3)(A6); n i is the number of patches in the landscape of patch type i (A2); E * is the total length (m) of edge in landscape (A4); e ik is the total length (m) of edge in landscape between patch types i and k, m is the number of patch types present in the landscape (A5); p ij is the perimeter of patch ij in terms of number of cell surfaces (A6).