Mapping and Characterizing Spatiotemporal Dynamics of Impervious Surfaces Using Landsat Images: A Case Study of Xuzhou, East China from 1995 to 2018

As an effective indicator of urbanization, impervious surfaces play a significant role in urban planning and ecological protection. It is, therefore, important to characterize impervious surfaces in urban geographical studies. As a key city in East China, Xuzhou has experienced rapid urbanization in recent decades and is now becoming an environmentally friendly city. To better understand the spatiotemporal heterogeneity of Xuzhou’s urban development, we extracted its impervious surfaces from Landsat images of 1995, 2003, 2010, and 2018 by a linear spectral mixture analysis. Then, a range of complementary methods including landscape indices, profile lines, median centers, standard deviational ellipses, and spatial autocorrelation were adopted to analyze the landscape pattern and expansion of impervious surfaces on both city and district scales. Results show that (1) there was a constant impervious surface expansion, originating in downtown Xuzhou; (2) promoting ecological protection in urban areas fragmented impervious surfaces with increasing heterogeneity and diversity overall; and (3) expansion directions and rates of impervious surfaces varied with district and town, and the central urban area expanded towards east and southeast, which could be related to their own resources and governmental policies. Findings from this study provide useful insights into urban planning of this economically prospective region.


Introduction
As population continues to increase, urbanization has become a universal and important phenomenon in the world. Urban expansion is one of the basic characteristics of urbanization, impacting urban ecosystem and changing urban land use structure and landscape pattern [1]. As such, appropriate methods are needed to accurately reflect the spatiotemporal variation of urban expansion, to understand the dynamics of regional urbanization.
A number of methods for mapping urban expansion exist, such as lognormal regression model [2], self-organizing map approach [3], landscape index [4,5], and (un-)supervised classification [6][7][8][9][10][11]. As the major component of land cover in urban areas, impervious surface (IS), i.e., artificial materials which water cannot penetrate, provides a key piece information on urban ecosystems and urban  Table 1), nine in Yunlong, ten in Jiawang, 11 in Gulou, 14 in Quanshan, and 24 in Tongshan district. The purple lines indicate the profile lines used in Section 3.3. The orange line in Tongshan district illustrates the division of this district used in Section 3.4.   Table 1), nine in Yunlong, ten in Jiawang, 11 in Gulou, 14 in Quanshan, and 24 in Tongshan district. The purple lines indicate the profile lines used in Section 3.3. The orange line in Tongshan district illustrates the division of this district used in Section 3.4.

Remote Sensing Images
Mapping impervious surface change of the study area provides a basis for understanding the pattern of urban development. Considering the rapid socioeconomic development of Xuzhou over the past decades, and the image quality and time interval of satellite remote sensing data, multi-temporal satellite image data from Landsat 5 TM (Thematic Mapper), Landsat 7 ETM+ (Enhanced Thematic Mapper Plus), and Landsat 8 OLI (Operational Land Imager), acquired in 1995, 2003, 2010, and 2018 (Table 2), were used to extract impervious surfaces and to identify the pattern, direction, and rate of impervious surface expansion.   Due to its large geographical size, the study area was imaged in two scenes (path 121/row 36, path 122/row 36). For each of the four years, it would be ideal if these two images were acquired at the same date and time. This is, however, difficult to achieve given image availability and image quality. As such, for each year, we collected two images whose acquisition data difference was within one week. This ensures that the difference was limited and would not have a significant impact on the results. In addition, as lush vegetation in summer and bare land in winter may interfere with the extraction of impervious surfaces, we considered images acquired in spring (March and April) more suitable for the study. All the remote sensing data were freely downloaded from the United States Geological Survey website (USGS, https://earthexplorer.usgs.gov/).
Prior to mapping impervious surfaces, preprocessing operations, such as atmospheric correction (FLAASH), geo-referencing (image to image), and seamless mosaic (in which a seamline network is formed by seamlines using area Voronoi diagrams with overlaps [43]) was applied to these images. Then, the study area was extracted from the mosaicked images using the vector data representing the extent of the study area.

Linear Spectral Mixture Analysis
The linear spectral mixture analysis (LSMA) was used to map impervious surface fractions from multi-temporal Landsat images. In this method, it is assumed that a pixel's spectral signal is a linear combination of the endmember spectral signals of each component [18]. The fraction of each endmember can be obtained by applying the least square model according to their own spectral characteristics [18]. In a fully contained linear spectral mixture analysis, two conditions must be met: (1) The sum of the fractions of all endmembers in each pixel is 1, and (2) the fraction of each endmember is non-negative. This is given as follows [18,28]: where R b is the reflectance for each band b in a pixel; f i is the fraction of endmember i in a pixel; R i,b is the reflectance of endmember i in band b for that pixel; and e b is the residual. The fitness of the fully constrained LSMA is assessed by the root mean square (RMS), which is calculated from e b .
The formula for RMS is as follow [18]: where n is the number of bands used in the LSMA. Water bodies should first be removed from the images to avoid noise as it would reduce unmixing errors caused by low-albedo impervious surfaces [44]. Existing methods of extracting water include principal component analysis (PCA) [45], spectral indices [46], and (un-)supervised classification [44]. Due to its simplification and accuracy, the modified normalized difference water index (MNDWI) [46], a widely used spectral index, was adopted to obtain and mask water body from the remote sensing data. Then, we adopted brightness normalization to decrease intra-class variability [18,19,30], and the minimum noise fraction (MNF) was adopted to separate the noise in the data and improve the quality of the endmember selection [18]. Endmember selection is critical to the accuracy of the model, and usually, three to four endmembers are appropriate [47]. Based on the VIS model [31], four endmembers were selected from the images by trial and error, namely vegetation, bare soil, high-and low-albedo impervious surfaces [18]. For each year, four fraction maps and one RMS map were produced through the LSMA, among which the high-and low-albedo impervious surface fraction maps were added together to generate an impervious surface fraction map.
The LSMA fitness was first assessed by RMS. If RMS is within 0.02, the LSMA results would be considered satisfactory [18]. To further assess the accuracy of impervious surface mapping, it is necessary to investigate the agreement between the mapped and real impervious surface fractions. Although the confusion matrix and Kappa coefficient are the commonly used accuracy assessment approaches [48,49], they are suitable for categorical variables, such as land cover classes from image classification. Impervious surface fractions are, however, continuous variables, and it would be more interesting to compare the mapped impervious surface fractions with real ones which were obtained from image digitalization and evaluate how close they would be [50,51]. In this study, a total of 80 validation samples, each sample containing 3 × 3 pixels (i.e., 90 m × 90 m), were randomly generated on the non-water area on the 2018 Landsat 8 image of the study area. Then, all the validation samples were overlaid on the high-resolution satellite images of 2018 in Google Earth Pro. The area of the digitized impervious surface in a given sample was divided by the area of the sample (i.e., 8100 m 2 ) as the true impervious surface proportion in the sample. The linear goodness of fit (R 2 ) between the impervious surface fraction in the LSMA and the true impervious surface proportion and the RMSE (root mean square error) for the 2018 LSMA were computed and assessed. It would be ideal if we could perform independent accuracy assessments for the modelled impervious surface fractions of each period. However, constrained by the poor image quality of historical Google Earth images and consequently, difficulty in digitalizing or visually interpreting them, we did not assess the extracted impervious surface fractions of 1995, 2003, and 2010. However, as all the images were Landsat data and processed with the same procedures, we assumed that similar accuracies were obtained for the LSMA results of 1995, 2003, and 2010. A similar assumption was also made by Cui et al. (2018) when they assessed their multi-temporal image classification results obtained from a supervised classifier [52].

Landscape Indices
Landscape indices are important indicators for quantitatively characterizing the composition and spatial distribution of the landscape pattern, which can reveal the structural changes in land use practices [53,54]. They might be flawed by spatial resolutions; thus, a proper selection of landscape indices is important for landscape studies [55]. Although there are some new measures of landscape fragmentation providing a good free-of-pixel size solution [56], we decided here to adopt more commonly used landscape indices for a detailed analysis (Table 3). They are able to detect changes at both class and landscape levels, not only the density, shape, and area of the patches (Patch density (PD), landscape shape index (LSI), largest patch index (LPI)) but their contagion, connectivity, and diversity (Aggregation index (AI), patch cohesion index (COHESION), Shannon's diversity index (SHDI), Shannon's evenness index (SHEI)). Diversity metrics (SHDI, SHEI) were used to reflect the diversity and heterogeneity of the landscape. The majority of these indices were also used in previous similar studies, such as Sha and Tian (2010) [57]. Note that these indices were calculated in Fragstats 4.2 based on the grading of impervious surface fractions (Section 4.1). Table 3. Landscape indices used in this study [58]. It is expressed such that an even distribution of area among patch types results in maximum evenness. Evenness is the complement of dominance

Profile Line Analysis
Profile line analysis was performed to analyze the directional and spatial variation of impervious surface fractions along these lines. As Pengcheng Square is in downtown Xuzhou as well as one of the most famous landmarks of Xuzhou, four profile lines in different directions-west-east, north-south, northwest-southeast, and southwest-northeast (Figure 1b), all through Pengcheng Square, were, therefore, adopted to measure the spatial difference of the impervious surface in different directions.

Median Center and Standard Deviational Ellipse
To further measure the direction of urban expansion, the median center and standard deviational ellipse (SDE) methods were used in this study. The median center is a measure of the central tendency, identifying the location where all other features are the least distanced in the data set [14]. The calculated median centers and their tracks on different scales reveal in which direction impervious surfaces expand.
The standard deviational ellipse is used to summarize spatial characteristics, such as central trends, dispersion, and directional trends, of all geographical features [39]. The semi-long axis and the semi-short axis of the ellipse indicate the direction and range of data distribution, respectively. Their ratio reflects the clustering or dispersion of impervious surfaces in a particular spatial direction. The greater the ratio is the more obvious the directionality of the data is, and vice versa. A long-short axis ratio equaling one indicates no directional characteristics of the data [39]. The rotation angle is the measure of the rotation of the semi-long axis, reflecting the change direction of impervious surfaces. If the features are densely centered and less dense toward the periphery, the SDE would encompass approximately 68% of the features. Identifying the positional variation of the SDEs of impervious surfaces allows one to see if the spatial distribution of impervious surfaces is elongated and shows a particular direction of impervious surface expansion, which can enhance the understanding of spatiotemporal characteristics of urban development.
Based on the impervious surface fraction maps, impervious surface binary maps (i.e., classifying the impervious surface fraction maps into impervious and pervious surfaces) were produced following the method described by Jia et al. [59]. Then, the median centers and standard deviational ellipses on different scales were calculated for each year.
It is worth noting that Tongshan district would be divided into two parts for median center and standard deviational ellipse analyses as it is by far larger than the other districts (Section 4.2.3). Based on the border of the CUA (Section 2.1), Tongshan district was divided into Tongshan North and Tongshan South (see the division line in Figure 1b) in this study for computing the median centers and SDEs of impervious surfaces.

Spatial Autocorrelation Analyses
Spatial autocorrelation refers to the interdependence between the attribute values of spatial units and the same attribute values of adjacent spatial units. It includes global and local spatial autocorrelation. Global spatial autocorrelation describes the average degree of association, spatial distribution patterns, and their significance among all geographic units in the study area, while local spatial autocorrelation can identify the aggregation and differentiation characteristics of local spatial features [52]. Moran's I and local Moran's I were used to describe the global and local spatial autocorrelation, given by [52]: where n is the total number of geographic units in the study area; x i and x j are the geographic units i and j, respectively (where i = j); x is the average value of all the geographic units; and w ij is the spatial weight matrix assigned to the pair of geographic units i and j (when i and j are topologically adjacent with a common point or a common edge, w ij is defined as 1, otherwise w ij is 0). A positive Moran's I value indicates a positive spatial autocorrelation (shown as spatial clusters, i.e., high-high clustering or low-low clustering), while a negative one indicates a negative spatial autocorrelation (shown as spatial outliers, i.e., the juxtaposition of high values next to low values). A zero value means that the variable obeys a stochastic spatial distribution [52]. The larger the absolute value, the stronger the spatial autocorrelation, and vice versa. Based on the Moran's I values, all geographic units are shown in a Moran's I scatter plot which has four quadrants [52]: high-high (HH), low-low (LL), high-low (HL), and low-high (LH). Corresponding to the Moran's I scatter plot, four types of geographic units, high-high, low-low, high-low, and low-high, are included in the LISA (local indicators of spatial association) cluster map, which only shows the geographic units that passed the significance test.
In this study, the spatial autocorrelation of impervious surface expansion in the study area was measured by the expansion speed index (ESI) to explore its spatial differences. Based on the impervious surface binary maps, ESI is defined as the ratio of the increased area of impervious surfaces to the original total impervious surface area within a given study area during a study period. It is given as follows [60]: where ∆S ij is the increased area of impervious surfaces in a geographical unit from period i to period j; S i is the total area of impervious surfaces in the geographical unit during period i; and ∆t is the period of time.

Impervious Surface Mapping
The average RMS over the images of LSMA was 0.017 for 1995, 0.019 for 2003, 0.020 for 2010, and 0.018 for 2018, all within 0.02. The accuracy assessment based on Google Earth samples for the 2018 LSMA result shows the linear goodness of fit (R 2 ) was 0.86 and the RMSE was 0.115. As such, we considered our LSMA results for 2018 as well as 1995, 2003, and 2010 were reliable and could be used for further analyses.
The impervious surface fraction maps obtained from the LSMA for different years are shown in Figure 2. It is clear that impervious surface fractions increased significantly from 1995 to 2018. Growing impervious surface fractions were observed in the downtown area and later, also in the area around it. However, impervious surface fractions of Tongshan district in the southeast and Jiawang district in the northeast remained relatively low.
where ∆ is the increased area of impervious surfaces in a geographical unit from period to period ; is the total area of impervious surfaces in the geographical unit during period ; and ∆ is the period of time.

Impervious Surface Mapping
The average RMS over the images of LSMA was 0.017 for 1995, 0.019 for 2003, 0.020 for 2010, and 0.018 for 2018, all within 0.02. The accuracy assessment based on Google Earth samples for the 2018 LSMA result shows the linear goodness of fit ( ) was 0.86 and the RMSE was 0.115. As such, we considered our LSMA results for 2018 as well as 1995, 2003, and 2010 were reliable and could be used for further analyses.
The impervious surface fraction maps obtained from the LSMA for different years are shown in Figure 2. It is clear that impervious surface fractions increased significantly from 1995 to 2018. Growing impervious surface fractions were observed in the downtown area and later, also in the area around it. However, impervious surface fractions of Tongshan district in the southeast and Jiawang district in the northeast remained relatively low.
Although used to be a mining city, Xuzhou has successfully transformed from rapid urbanization to environmental friendliness, which could be revealed by its vegetation cover changes. The mean vegetation fraction was 35.05% for 1995, 27.76% for 2003, 29.49% for 2010, and 43.04% for 2018 respectively, showing a significant rise after an early decline.  Although used to be a mining city, Xuzhou has successfully transformed from rapid urbanization to environmental friendliness, which could be revealed by its vegetation cover changes. The mean vegetation fraction was 35.05% for 1995, 27.76% for 2003, 29.49% for 2010, and 43.04% for 2018 respectively, showing a significant rise after an early decline.
To further analyze the spatiotemporal pattern of urban impervious surfaces, impervious surface fractions were divided into five levels: 0-0.2 for low-density impervious surfaces, 0.2-0.4 for medium-low-density, 0.4-0.6 for medium-density, 0.6-0.8 for medium-high-density, and 0.8-1 for high-density. The levels of the impervious surface fractions for the four years are shown in Figure 3. The number of pixels decreased remarkably for the low-density impervious surfaces over time but increased gradually for the other levels. To further analyze the spatiotemporal pattern of urban impervious surfaces, impervious surface fractions were divided into five levels: 0-0.2 for low-density impervious surfaces, 0.2-0.4 for medium-low-density, 0.4-0.6 for medium-density, 0.6-0.8 for medium-high-density, and 0.8-1 for high-density. The levels of the impervious surface fractions for the four years are shown in Figure 3. The number of pixels decreased remarkably for the low-density impervious surfaces over time but increased gradually for the other levels.

Based on Landscape Indices
Based on the five levels of impervious surface fractions (Section 4.1), landscape indices characterizing impervious surfaces were calculated in Fragstats 4.2 using the method described in Section 3.2. Results are shown in Figure 4.

Based on Landscape Indices
Based on the five levels of impervious surface fractions (Section 4.1), landscape indices characterizing impervious surfaces were calculated in Fragstats 4.2 using the method described in Section 3.2. Results are shown in Figure 4.
The PD values gradually increased for all the IS fraction levels except for that for the medium-low-density level (0.2-0.4) which remained the highest with the medium-density level (0.4-0.6) levels (Figure 4a). The LSI values also showed a rising trend for all IS fraction levels except for the low-density level (0-0.2) which saw a gentle decline in 2018 (Figure 4b). Both AI and COHESION values showed similar changes over time (Figure 4c,d). The values for the low-density (0-0.2) and medium-low-density (0.2-0.4) levels decreased gradually while those for the other IS fraction levels rose before going down. The LPI values for the low-density level (0-0.2) decreased remarkably while little change was observed for the other levels ( Figure 4e). Regarding the values of landscape-level SHDI and SHEI, both of them increased steadily over the periods (Figure 4f). The PD values gradually increased for all the IS fraction levels except for that for the mediumlow-density level (0.2-0.4) which remained the highest with the medium-density level (0.4-0.6) levels ( Figure 4a). The LSI values also showed a rising trend for all IS fraction levels except for the lowdensity level (0-0.2) which saw a gentle decline in 2018 (Figure 4b). Both AI and COHESION values showed similar changes over time (Figure 4c and 4d). The values for the low-density (0-0.2) and medium-low-density (0.2-0.4) levels decreased gradually while those for the other IS fraction levels rose before going down. The LPI values for the low-density level (0-0.2) decreased remarkably while little change was observed for the other levels ( Figure 4e). Regarding the values of landscape-level SHDI and SHEI, both of them increased steadily over the periods (Figure 4f).

Based on Profile Line Analysis
Analyses in Sections 4.1 show that impervious surface expansion varies within the study area. This might be further revealed through profile lines analysis. Along the four profile lines through Pengcheng Square (Figure 1b

Based on Profile Line Analysis
Analyses in Section 4.1 show that impervious surface expansion varies within the study area. This might be further revealed through profile lines analysis. Along the four profile lines through Pengcheng Square (Figure 1b), we extracted impervious surface fractions for each year ( Figure 5).
Along the W-E profile line traveling from Dapeng town and ending with Daxu town in Tongshan district through Quanshan, Gulou, Yunlong, and Gulou district (Figure 5a), impervious surface fractions were lower in the east than those in west. The N-S profile line also started and ended in the Tongshan district through the districts of Gulou, Yunlong, and Quanshan (Figure 5b) reveals that the impervious surface fractions of the CUA were higher than the other areas along the line. The NW-SE profile line that travels from Huangji town to Fangcun town (Figure 5c) also shows the highest impervious surface fractions were focused in Yunlong and Gulou districts, and the lowest in the southeast of Tongshan district. Regarding the SW-NE profile line going through Hanwang town of Tongshan district, Quanshan, Gulou, Tongshan, and Jiawang districts (Figure 5d), impervious surface fractions were clearly higher the CUA than the other areas. From a spatial perspective, the CUAs were constantly observed with the highest impervious surface fractions; from a temporal perspective, impervious surface fractions increased in general over the 23 years, though differently.

Based on Median Centers and Standard Deviational Ellipses
Following the method described in Section 3.4, the median centers of impervious surfaces were computed, and their tracks were illustrated for both the entire study area ( Figure 6) and individual districts (Figure 7). On the city scale, the median centers of the impervious surfaces were located at  (Table 4), and shifted nearly horizontally over time. During the 2003 to 2010 period, the annual average movement to the east accelerated. However, the overall moving distance from 1995 to 2018 was relatively small for the entire area ( Figure 6), and the overall expansion direction is not significant. On the district scale, however, the median centers of the impervious surfaces in Gulou, Yunlong, and Quanshan districts showed obvious expansion directions, towards the east (Figure 7b), southeast (Figure 7c), and northwest ( Figure 7f). Tongshan district shows internal differences, with the expansion directions and distances in Tongshan North were different from those in Tongshan South (Figure 7, Table 4).
perspective, impervious surface fractions increased in general over the 23 years, though differently.

Based on Median Centers and Standard Deviational Ellipses
Following the method described in Section 3.4, the median centers of impervious surfaces were computed, and their tracks were illustrated for both the entire study area ( Figure 6) and individual districts (Figure 7). On the city scale, the median centers of the impervious surfaces were located at 117.264° E-117.309° E, 34.311° N-34.316° N during the 23-year period (Table 4), and shifted nearly horizontally over time. During the 2003 to 2010 period, the annual average movement to the east accelerated. However, the overall moving distance from 1995 to 2018 was relatively small for the entire area (Figure 6), and the overall expansion direction is not significant. On the district scale, however, the median centers of the impervious surfaces in Gulou, Yunlong, and Quanshan districts showed obvious expansion directions, towards the east (Figure 7b), southeast (Figure 7c), and northwest (Figure 7f). Tongshan district shows internal differences, with the expansion directions and distances in Tongshan North were different from those in Tongshan South (Figure 7, Table 4).     The standard deviational ellipses (SDEs) are shown in Figure 8. For the entire study area, the rotation angles ranged from 96.539 • to 109.964 • , and the long-short axis ratios were all close to one during the 1995 to 2018 period ( Table 5), indicating that the impervious surfaces had no clear expansion direction. The SDEs on the district scale were different from those on the city scale ( Figure 8). The long-short axis ratio of each district was significantly larger than that of the entire study area (Table 5), reflecting the obvious expansion direction on the district scale. As the major components of the CUA, Gulou, Yunlong, and Quanshan districts presented expansion towards the east (Figure 8b), southeast (Figure 8c), and northwest ( Figure 8f). Such observations on both the city and district scales were consistent with those of the median centers in Figure 7 and Table 4.

of 23
The standard deviational ellipses (SDEs) are shown in Figure 8. For the entire study area, the rotation angles ranged from 96.539° to 109.964°, and the long-short axis ratios were all close to one during the 1995 to 2018 period ( Table 5), indicating that the impervious surfaces had no clear expansion direction. The SDEs on the district scale were different from those on the city scale ( Figure  8). The long-short axis ratio of each district was significantly larger than that of the entire study area (Table 5), reflecting the obvious expansion direction on the district scale. As the major components of the CUA, Gulou, Yunlong, and Quanshan districts presented expansion towards the east ( Figure  8b), southeast (Figure 8c), and northwest ( Figure 8f). Such observations on both the city and district scales were consistent with those of the median centers in Figure 7 and Table 4.      (Figure 9), showing that there were positive autocorrelations in the expansion rates of impervious surfaces. In addition, the Moran's I gradually increased over time (Figure 9). The towns of ESIs in the three periods were mainly distributed in the HH and LL quadrants, indicating that spatial homogeneity and aggregation pattern of expansion rate were significant. In addition, the towns in the LL quadrants is aggregate, while the towns in the HH quadrants are scattered.
In addition, the LISA cluster maps ( Figure 10) were also generated to show the type of towns which passed the significance test. We observed that the HH and LL clustering dominated in the three periods. In the 1995 to 2003 period (Figure 10a), the HH and LH clustering occurred around Pengcheng Square and the LL clustering in Jiawang and Tongshan districts. However, there seemed an opposite trend that LL and HL clustering are distributed around Pengcheng Square for the 2003 to 2010 period (Figure 10b). In the 2010 to 2018 period (Figure 10c), both the HH and LL clustering expanded compared to the previous period.
both global and local spatial autocorrelation analyses were performed to explore their spatial correlations.
The Moran's I of ESIs was 0.2748, 0.4461, and 0.6748 in [1995][1996][1997][1998][1999][2000][2001][2002][2003][2003][2004][2005][2006][2007][2008][2009][2010], and 2010-2015, respectively (Figure 9), showing that there were positive autocorrelations in the expansion rates of impervious surfaces. In addition, the Moran's I gradually increased over time (Figure 9). The towns of ESIs in the three periods were mainly distributed in the HH and LL quadrants, indicating that spatial homogeneity and aggregation pattern of expansion rate were significant. In addition, the towns in the LL quadrants is aggregate, while the towns in the HH quadrants are scattered.  In addition, the LISA cluster maps ( Figure 10) were also generated to show the type of towns which passed the significance test. We observed that the HH and LL clustering dominated in the three periods. In the 1995 to 2003 period (Figure 10a), the HH and LH clustering occurred around Pengcheng Square and the LL clustering in Jiawang and Tongshan districts. However, there seemed an opposite trend that LL and HL clustering are distributed around Pengcheng Square for the 2003 to 2010 period (Figure 10b). In the 2010 to 2018 period (Figure 10c), both the HH and LL clustering expanded compared to the previous period.

Interpretation and Discussion
In this study, we first mapped urban impervious surfaces of Xuzhou by the LSMA from multitemporal Landsat image data of 1995, 2003, 2010, and 2018. As urban expansion is a complex process involving landscape pattern, expansion direction, and rate difference, various methods including landscape indices, profile lines, median centers, standard deviational ellipses, and spatial autocorrelation were adopted to characterize the above aspects and comprehensively reflect urban dynamics. The interpretation of the results and their implications are detailed below.

Interpretation and Discussion
In this study, we first mapped urban impervious surfaces of Xuzhou by the LSMA from multi-temporal Landsat image data of 1995, 2003, 2010, and 2018. As urban expansion is a complex process involving landscape pattern, expansion direction, and rate difference, various methods including landscape indices, profile lines, median centers, standard deviational ellipses, and spatial autocorrelation were adopted to characterize the above aspects and comprehensively reflect urban dynamics. The interpretation of the results and their implications are detailed below.

Overall Dynamics
Impervious surfaces of Xuzhou rose dramatically with rapid urbanization in the recent two decades, as there exists a strong relationship between impervious surface fractions and urban development [61]. The downtown area has always had relatively high impervious surface fractions due to its high socioeconomic level and building density (Figure 2). The rapid development of the CUA has been driving the development of the whole city over time [62] and shows a positive impact on its neighboring towns (Figure 2), a pattern observed in other cities in Asia, such as Tokyo [63] and Shanghai [64]. The transformation of impervious surfaces from low density to high density is a typical feature of rapid urbanization ( Figure 3) [65]. In addition, the development of the new urban area (SE Yunlong district) and the construction of transportation facilities in Xuzhou, such as the subway, further promote the overall urbanization, leading to increased impervious surface fractions and a complicated and fragmented trend in the spatial pattern [54], showing the agglomeration effect of urban impervious surfaces. This is because of the implementation of urban renewal projects in Xuzhou-scattered and chaotic buildings have been mostly replaced with continuous residential, commercial, and industrial facilities, leading to increased impervious surface fractions. As there is a negative correlation between impervious surfaces and vegetation [31], low impervious surface fractions were observed in the southeast of the study area where Lvliang Hill Scenic spots was located with abundant vegetation. Balanced development of impervious surface expansion and ecological environment protection contributes to an increased landscape diversity.
Along with urban expansion, vegetation has also flourished in recent years. In the early period, vegetated land was converted to built-up land for constructing residential and commercial buildings. As ecological protection consensus started to emerge, efforts, such as afforestation and ecological restoration for coal mining subsidence areas, have brought more greenness back to Xuzhou in recent years. These measures were also proven effective in other cities [66].

Landscape Pattern
An increase in PD indicates increasing landscape heterogeneity of impervious surfaces (Figure 4a). The advancement of urbanization has made impervious surfaces expand outside with a large number of residential areas scattered in the suburbs, leading to an increasing fragmentation trend of impervious surfaces. The increase in the LSI (Figure 4b) also indicates that the spatial pattern of impervious surfaces is irregular and complicated. In addition, the establishment of industrial zones, the development of the new urban areas (SE Yunlong district) and the construction of transportation facilities in Xuzhou, such as subways, further promoted the overall urbanization, making the original landscape pattern more complex and fragmented [54,65]. The changes in AI and COHESION (Figure 4c,d) indicate that the low-and medium-low-density impervious surfaces exhibited a dispersive expansion and the connectivity was reduced. This is because the promotion of eco-environmental protection, such as the construction of scenic spots, made the impervious surfaces with low density become fragmented and irregular. While the impervious surfaces with high density showed a filled expansion and enhanced connectivity, which is closely related to the transformation of the old city in Xuzhou, the LPI of the low-density impervious surface dramatically reduced (Figure 4e), indicating that urbanization made its dominance significantly reduce, and part of it was gradually occupied by high-density [65], which is the inevitable result of urban expansion.
At the landscape level, gradually increased SHDI and SHEI (Figure 4f), suggest enhanced heterogeneity and landscape diversity. This is because urbanization has reduced the dominance of agricultural land and vegetation by converting part of them into impervious surfaces. Meanwhile, there are some open green spaces in the study area, such as Yunlong Mountain scenic area, Quanshan forest park, and Jiulihu wetland park, and large-area water landscapes, such as Yunlong Lake and Dalong Lake. These permeable landscapes are important for maintaining the original components and the quality of urban ecological environment [67,68]. Coordinated development of impervious surface expansion and eco-environment protection means the landscape types are in a relatively balanced state and increases landscape diversity.

Expansion Direction
Profile lines, median centers, and SDEs are widely used to characterize urban expansion [14,39]. While Xu et al. (2018) [39] found that the expansion direction of Guangzhou was significantly contrasting on the whole and local region scale, we found Xuzhou featured different expansion directions than Guangzhou.
In this study, the impervious surfaces did not exhibit an obvious expansion direction on the city scale ( Figures 6 and 8). The annual average rate of expansion to the east has increased during the 2003 to 2010 period (Table 4) probably resulted from a governmental policy-the Xuzhou Eleventh Five-Year Plan (2006-2010) had clearly stated Xuzhou would expand eastward. However, this was not significant for the entire study area from 1995 to 2018. On the district scale ( Figures 5-8), Quanshan, Gulou, and Yunlong districts, the major components of the CUA, have increasingly high impervious surface fractions, which is related to their urbanization that started early and developed rapidly. The impervious surfaces of Quanshan district expanded towards the northwest before the Pangzhuang mine was closed in Pangzhuang town in the NW Quanshan district. While coal mining related subsidence was a serious issue in Pangzhuang town, the government encouraged the expansion of the CUA towards east and southeast, e.g., building a key high-speed railway station in east Gulou district and creating the Xuzhou national economic development zone covering east CUA, relocating government offices to the new urban area in SE Yunlong district. These expansion directions are consistent with the requirements of the Xuzhou General Urban Planning (2007-2020) that built-up land in the CUA should be developed in the east and southeast, be controlled in the northwest, and strictly restricted in the southwest [69].

Expansion Rate
The expansion rates of impervious surfaces during the two decades were shown to have significant spatial autocorrelation, and the spatial agglomeration gradually increased with the advancement of urbanization (Figure 9). HH and LL clustering of impervious surface expansion dominated in Xuzhou, the central city of the Huaihai Economic Zone (Figure 9), as urban development in economic regions tends to be homogeneous [70].
The urbanization of the downtown area started early, showing HH and LH clustering around Pengcheng Square in the 1995 to 2003 period (Figure 10a). Note that the clustering only implies a relative change, rather than an absolute value difference [14]. LH clustering does not necessarily mean that the towns in this type have a low expansion rate, but tells us that their neighboring towns' impervious surfaces expand faster. In the periods of 2003 to 2010 and 2010 to 2018, opposite spatial distributions were observed with the LL and LH clustering occurring in the downtown area (Figure 10b,c), which is due to the fact that the development of downtown area slowed down relatively after approaching saturation. Likewise, the town shown as HL clustering expanded relatively faster because of their proximity to the saturated downtown area (Figure 10b). In addition, Tongshan and Jiawang, which are not included in the CUA, developed late, but their urbanization has accelerated in recent years (Figure 10b,c). Tongshan district developed rapidly, upgrading from a rural county to an urban district. This explains why LL clustering appeared in the downtown area, and HH clustering were in the suburbs. However, because of unevenly natural resource distribution [71] and the effect of government policies, expansion rate varied in the towns largely within Tongshan district (Table 4), with the HH clustering in different locations (Figure 10c).

Innovation and Limitation
Studies of the spatiotemporal dynamics of impervious surfaces on either region or city scale are important for understanding the overall pattern and development model [33][34][35][36]. However, these may fail to notice the variation among the areas of lower administrative levels in regional development [72][73][74]. A better characterization of the spatiotemporal dynamics of impervious surfaces entails perspectives on different scales in the case study of Xuzhou. In addition, although a single method may unravel a certain aspect of urban development [39,41], this is overshadowed by the complexity of urban development [75]. The combination of many different methods together with analyses on different scales, as demonstrated in our case study, provides a new perspective for comprehensive quantification of urban expansion.
However, we here would like to raise the issue of discriminating endmembers in the process of endmember selection in this study. While vegetation is most detectable and low-albedo impervious surfaces exist in large amounts in urban areas, it is challenging to identify high-albedo impervious surfaces, and bare soil endmembers as both of them are usually uncommon in cities. Moreover, because they are quite similar in spectral characteristics, high-albedo impervious surfaces could be easily confused with bare soil. As such, an increased resolution of image data and a more advanced endmember extraction procedure would improve the accuracy of mapping impervious surfaces.

Conclusions
This study characterizes and analyzes the spatiotemporal dynamics of impervious surfaces of Xuzhou from 1995 to 2018 on different scales using various and complementary methods. The key findings and main conclusions are summarized as follows: • Impervious surfaces increased obviously in the context of rapid urbanization, which changed urban landscape patterns. Impervious surfaces with high fractions were mainly concentrated in the downtown area, showing an expansion starting from the downtown area. Impervious surfaces were generally fragmented and irregular. Meanwhile, vegetation also flourished in recent years. Scientific urban planning promotes the balanced development of impervious surface expansion and ecological environmental protection, increasing the diversity of landscape. • Significant differences in the expansion direction of impervious surfaces existed in the entire study area and each district. The expansion direction of the study area was not obvious, while the districts within the CUA shows clear expansion directions towards the east and southeast, which is consistent with the general urban planning. Therefore, more importance should be placed on the urban planning and policy guidance to stimulate and regulate the overall orderly urban development.

•
Expansion rates of impervious surfaces showed a significant spatial agglomeration, which increased gradually and varied with the town. The urbanization of the downtown area started early and has gradually become saturated, while the non-CUA accelerated its development with the large internal differences. This suggests that resource distribution and government policies affect urban expansion rates.
Findings from this study will be helpful to understand the tempo-spatial variation in urban expansion, providing a basis for analyzing its environmental effects and exploring its driving mechanism. Moreover, this will contribute to decision-makers in formulating urban planning policies and provides a reference for similar studies.