Monitoring Three-Decade Expansion of China’s Major Cities Based on Satellite Remote Sensing Images

As the largest developing country, China has experienced dramatic urban expansion since the “reform and opening-up” policy started at the end of the 1970s. In this paper, we monitor three decades of urban expansion in China’s 36 major cities, based on the spectral mixture analysis of remotely sensed satellite images. The results demonstrated that these major cities have expanded by 5.85 times from 1986 to 2015, with 15.51 km2 average expansion area per city per year. We found the urban expansion trajectories showed three different modes, i.e., exponential, linear and s-shaped, which were closely related to the city development level. In the old city zones, however, there was an interesting common tendency of the impervious surface area (ISA) first increasing and then decreasing, which could be largely attributed to the phenomenon of urban village reconstruction in China’s cities. Based on the Granger Causality Test (GCT), the interaction between urban ISA and Gross Domestic Product (GDP) per capita (GDPPC) suggested that the former was the driver of the latter. Meanwhile, taking the Yangtze River as the division between north and south China, there exists a north–south territorial differentiation for the interaction between ISA and total population at the year-end (TP).


Introduction
Urbanization can be expressed as a shift of population from rural to urban areas, predominantly resulting in a transformation from rural land to urban land [1][2][3][4]. According to the United Nations, by 2014, 54% of the world population had lived in an urban area [5], and this is predicted to reach 66% by 2050. The growth of the urban population inevitably results in urban area expansion. It has been reported that the global urban area is likely to increase by 1.2 million km 2 by 2030 [6]. Although urban areas represent a small fraction of the Earth's surface, they have global impacts, including natural land loss, food waste, ecological degradation, environmental pollution, and so on [4,[6][7][8][9][10].
Remotely sensed data, with superiority of short update cycle, continuous functioning and large coverage [11], has become an important means for urban studies around the world, especially for large-scale monitoring. At a global scale, the global urban areas have been measured at coarse spatial

Data Preprocessing
Radiometric calibration, geometric rectification, mosaicking and clipping of images were all conducted. In detail, for each city, if the images were not georeferenced with each other, the image in the first study year was used as the reference to perform image-to-image registration, and then every image was clipped by the administrative district boundary. Image mosaicking was applied in Changchun and Shanghai. Moreover, to reduce the influence of atmospheric effects, the Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) atmospheric correction algorithm was applied to all experimental images.
To acquire more experimental data, an image recovery [61,62] method was applied in this research. For the study period of 2003-2015, if the classification results for the available cloud-free Landsat images were not acceptable, and the time span between two adjacent study years was greater than three years, we filled the gaps in the SLC-off imagesfollowing the image recovery method in [63,64] (software could be downloaded from The data used included two main parts. (1) Satellite data. Cloud-free Landsat images (digital number data) covering the 36 major cities in mainland China were obtained from the United States Geological Survey (USGS) website (http://glovis.usgs.gov/). One image was selected for one year. From Table S1, it could be found that more than 17 images were collected for the past three decades for each city, with a time span between two adjacent study years less than three years. (2) Statistical data. The statistical data (GDP and total population at the year-end) were obtained from the China City Statistical Yearbook . The result of GDP divided by total population was regarded as GDP per capita.

Data Preprocessing
Radiometric calibration, geometric rectification, mosaicking and clipping of images were all conducted. In detail, for each city, if the images were not georeferenced with each other, the image in the first study year was used as the reference to perform image-to-image registration, and then every image was clipped by the administrative district boundary. Image mosaicking was applied in Changchun and Shanghai. Moreover, to reduce the influence of atmospheric effects, the Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) atmospheric correction algorithm was applied to all experimental images.
To acquire more experimental data, an image recovery [61,62] method was applied in this research. For the study period of 2003-2015, if the classification results for the available cloud-free Landsat images were not acceptable, and the time span between two adjacent study years was greater than three years,

Urban Area Extraction
Various definitions for urban area will result in diverse estimates [41]. Existing definitions for remote sensing data can be classified as two types. One is defined as all the built-up area within the administrative division [12,41,59,67]. The other is defined as the urban core area and its suburbs [39,44,58]. In this paper, we chose the latter definition of urban extent, which could be amply described as the administrative division, the area within the core built-up land with high-density ISA and those built-up clusters closely related to the urban core and dispersed in suburbs, including all impervious and pervious surfaces. The suburb was defined as a 0-5km buffer outside the urban core area [68,69]. Specifically, the new districts set up by government were considered as an urban area even though some of them were far away from the core built-up area. The method used to obtain the precise urban extent involved two steps.

Impervious Surface Retrieval
The urban area in this research was generated based on impervious surface. To acquire a precise ISA, LSMA [65,[70][71][72] was used in this manuscript. In previous research, it has been suggested that the ISA should be classified by their essential impervious materials [73]. Considering the ground truth of China, the ISA was classified as fiber cement board (roofs), cement (roads, bridges), color coated steel sheets (sheds), brick (roofs), and asphalt (roads). For each year, we collected the endmembers from the Landsat image corresponding to the year. More details about the endmembers selection can been found in Supplementary Information Text S1. With the extracted According to previous research, water bodies may be easily confused with the low-reflectance impervious surface [65]. The Normalized difference water index (NDWI) [66] was therefore used to exclude water bodies before conducting the urban area extraction. For each image, the NDWI threshold was set as 0,0.1,0.2,0.3,0.4, and 0.5, respectively. The best result with the least water bodies and the most impervious surface would be selected visually.

Urban Area Extraction
Various definitions for urban area will result in diverse estimates [41]. Existing definitions for remote sensing data can be classified as two types. One is defined as all the built-up area within the administrative division [12,41,59,67]. The other is defined as the urban core area and its suburbs [39,44,58]. In this paper, we chose the latter definition of urban extent, which could be amply described as the administrative division, the area within the core built-up land with high-density ISA and those built-up clusters closely related to the urban core and dispersed in suburbs, including all impervious and pervious surfaces. The suburb was defined as a 0-5 km buffer outside the urban core area [68,69]. Specifically, the new districts set up by government were considered as an urban area even though some of them were far away from the core built-up area. The method used to obtain the precise urban extent involved two steps.

Impervious Surface Retrieval
The urban area in this research was generated based on impervious surface. To acquire a precise ISA, LSMA [65,[70][71][72] was used in this manuscript. In previous research, it has been suggested that the ISA should be classified by their essential impervious materials [73]. Considering the ground truth of China, the ISA was classified as fiber cement board (roofs), cement (roads, bridges), color coated steel sheets (sheds), brick (roofs), and asphalt (roads). For each year, we collected the endmembers from Remote Sens. 2020, 12, 491 5 of 19 the Landsat image corresponding to the year. More details about the endmembers selection can been found in Supplementary Information Text S1. With the extracted endmembers, a fully constrained linear spectral unmixing method was applied to quantify the land use abundance [65,[70][71][72].
where R b is the reflectance for each band b in an image pixel; R i,b represents the reflectance of endmember i in band b for that pixel; f i stands for the fraction of endmember i; N is the numbers of endmembers; and e b is the residual. The final ISA image was obtained after adding together all the impervious material fraction images.
For each year, the generated ISA image would be compared with the original true-color satellite data, so that some other land covers confused with ISA would be excluded manually. In particular, we focused more on the "boundary" around urban and rural. Furthermore, for those years (mostly after 2000) we were able to find high-resolution images from Google Earth, and the generated ISA image would also be compared with the high-resolution images. Using root-mean-square error (RMSE, Equation (3)) asthe measurement index for the ISA fraction, we validated the accuracy for ISA fractions.
where V i,t is the "true" ISA composition fraction digitized from Google Earth for samplei; V i is the estimated ISA composition fraction for sample i; and N is the sample numbers.74 samples were selected for accuracy assessment: 39 samples for Wuhan city (13 samples for 2000, 13 samples for 2004 and 13 samples for 2013, Figure 3) and 35 samples for the other cities (one sample for one city in 2015, Table S4). The sample size was set as 50 × 50 pixels. The estimated abundance was calculated as the mean ISA abundance within the 50 × 50 pixels inthe ISA image, resulting in the tested ISA abundances varying from 0.05 to 0.97 (Table S4). The true ISA within each sample was visually interpreted from the Google Earth imagery. By dividing the ISA by 2.25 km 2 (50 × 50 pixels, i.e., 1500 m × 1500 m), the result could be regarded as the true ISA abundance. More description about samples selection can be found in Supplementary Information Text S2. Table S4). The sample size was set as 50 × 50 pixels. The estimated abundance was calculated as the mean ISA abundance within the 50 × 50 pixels inthe ISA image, resulting in the tested ISA abundances varying from 0.05 to 0.97 (Table S4). The true ISA within each sample was visually interpreted from the Google Earth imagery. By dividing the ISA by 2.25 km 2 (50 × 50 pixels, i.e., 1500 m × 1500 m), the result could be regarded as the true ISA abundance. More description about samples selection can be found in Supplementary Information Text S2.

Generation of the Urban Area
Corresponding to the definition of urban area in this paper, the urban area could be generated based on the high-density ISA images. Referring to previous research [71,74], the ISA fractions greater than 0.5 would be assigned as high-density ISA. Firstly, aggregation was applied on the ISA image. For each target pixel, searching in its eight-neighborhood, the target pixel and its neighborhood pixels whose ISA fraction was greater than 0.5 (i.e., high-density) were assigned as the same value ( Figure 4), which was a similar approach to the City Clustering Algorithm (CCA) [74][75][76]. The search was undertaken until all the pixels were tested, after which an indexed map could be generated with a unique index for each aggregated area ( Figure 5b). Secondly, with the help of the high-resolution images, the aggregated areas corresponding to the urban core, built-up clusters in suburb and the new districts set up by the government (e.g., Tianjin Binhai New District) were extracted. For example, from satellite images, we were able to find that the urban core for Wuhan city was located alongside the Yangtze River. Thus, the indexes for aggregated areas corresponding to the urban core could easily be found by visual interpretation. With those indexes extracted, the urban area could be picked up from Figure 3b,c by attribute extraction. Finally, the boundary of the aggregated built-up area was generated ( Figure 5d). The areas within the boundary were assigned as the urban area, including all the impervious and pervious surfaces except water bodies.

Generation of the Urban Area
Corresponding to the definition of urban area in this paper, the urban area could be generated based on the high-density ISA images. Referring to previous research [71,74], the ISA fractions greater than 0.5 would be assigned as high-density ISA. Firstly, aggregation was applied on the ISA image. For each target pixel, searching in its eight-neighborhood, the target pixel and its neighborhood pixels whose ISA fraction was greater than 0.5 (i.e., high-density) were assigned as the same value (Figure 4), which was a similar approach to the City Clustering Algorithm (CCA) [74][75][76]. The search was undertaken until all the pixels were tested, after which an indexed map could be generated with a unique index for each aggregated area ( Figure 5b). Secondly, with the help of the high-resolution images, the aggregated areas corresponding to the urban core, built-up clusters in suburb and the new districts set up by the government (e.g., Tianjin Binhai New District) were extracted. For example, from satellite images, we were able to find that the urban core for Wuhan city was located alongside the Yangtze River. Thus, the indexes for aggregated areas corresponding to the urban core could easily be found by visual interpretation. With those indexes extracted, the urban area could be picked up from Figure 3b,c by attribute extraction. Finally, the boundary of the aggregated built-up area was generated ( Figure 5d). The areas within the boundary were assigned as the urban area, including all the impervious and pervious surfaces except water bodies. When a target pixel is selected (the left image), a search is carried out within its eight-neighborhood until all the pixels greater than 0.5 are found (the middle image). Next, these pixels are assigned as the same number N. If this is the first search, then N would be 1; for the second search, N would be 2; and so forth.
Since all the impervious and pervious surfaces (water bodies excluded) within the urban boundary were marked as urban area, it could be pointed out that the land cover changes within the urban boundary would not affect the total area of urban extent. Thus, considering the progress of urbanization, we insisted the hypothesis that the urban area in a particular year would not be less When a target pixel is selected (the left image), a search is carried out within its eight-neighborhood until all the pixels greater than 0.5 are found (the middle image). Next, these pixels are assigned as the same number N. If this is the first search, then N would be 1; for the second search, N would be 2; and so forth.
Since all the impervious and pervious surfaces (water bodies excluded) within the urban boundary were marked as urban area, it could be pointed out that the land cover changes within the urban boundary would not affect the total area of urban extent. Thus, considering the progress of urbanization, we insisted the hypothesis that the urban area in a particular year would not be less than that in the former year. In practice, the urban area in this year can be expressed as the sum of the new urban area generated in this year and the urban area of the last year.  Since all the impervious and pervious surfaces (water bodies excluded) within the urban boundary were marked as urban area, it could be pointed out that the land cover changes within the urban boundary would not affect the total area of urban extent. Thus, considering the progress of urbanization, we insisted the hypothesis that the urban area in a particular year would not be less than that in the former year. In practice, the urban area in this year can be expressed as the sum of the new urban area generated in this year and the urban area of the last year.

Accuracy of Impervious Surface Area
The analysis for the 36 major cities were generally based on the ISA extracted in this paper. Thus, the spatial and temporal accuracy for ISA was required to be accurate and stable. Based on the aforementioned method of accuracy assessment, we were pleased to find that the extracted ISA was a series with stable accuracy.
In detail, from Table S4 (right part), taking Wuhan as the test area, there existed an average error of 0.05-0.07 for different years (0.05 for 2000; 0.07 for 2004; and 0.06 for 2013, average for three years: 0.06). Thus, the temporal consistency for the ISA was tested. For the spatial accuracy, we were able to find that the difference between estimated ISA abundance and true ISA abundance was mostly less than 0.1 with an average RMSE of 0.07 for all abundances (0.05-0.97, value interval circa 0.28) across all cities. As the spatial and temporal error level for ISA abundances was relatively low and stable, the value and temporal trends for the following extracted ISA and urban area were also acceptable.

Urban Expansion Rates and Modes
To qualify the urban expansion area and rates, by analyzing the three decades of satellite images, our results showed that almost all the cities experienced a dramatic expansion in this period. The urban area ratio (UAR) was calculated as the urban area in 2015 divided by the urban area in 1986. For those cities without clear images in 1986 or 2015, extrapolation or interpolation (for those cities with clear images before 1986 or after 2015) was used to obtain the urban area in 1986 or 2015. The UAR for all 36 cities was 5.85, and the expansion area was 15.51 km 2 per city per year. In addition, many of the cities showed an average growth rate (AGR , Table 1) [36,41] of around 4-7%. However, large variations in expansion area and AGR exist across China. For instance, Shenzhen (Figure 6b), as the "window" of China's reform and opening-up, possesses the highest AGR (12.23%) and the largest UAR (31.88), living up to the title of "a city rising in a night". Meanwhile, the lowest level of urban expansionwas found in Beijing (Figure 6a), with an AGR of 3.23% and an UAR of 2.59. This may be attributed to two reasons. Firstly, as the capital of China, Beijing had already developed to a relatively high level by the 1980s, and thus the available space for continued urban expansion was limited. Secondly, the stringent urban land policies in Beijing played an important role in restraining the rapid expansion. As for the expansion area (Figure 7), Shanghai tops the list, with 1716.03km 2 taken up by urban expansion. For one thing, the original urban area in Shanghai was large. For another thing, the Pudong district, which developed from the late 20th century and is regarded as the second national new district after Shenzhen, promoted the urban expansion in Shanghai. Meanwhile, Lhasa, restricted by natural conditions, shows the smallest increase in urban extent, with only 96.42 km 2 taken up by urban expansion. Overall, China has experienced dramatic urban expansion; however, influenced by a development basis, government policies, and natural conditions, the urban expansion shows great heterogeneity across cities.
To further explore the expansion patterns for different cities, a map illustrating the urban expansion trajectories is provided in Figure 8. For the city areas generated in the different study years, a linear function (linear expansion), exponential function (exponential expansion), and logistic function (s-shaped expansion) were applied to fit the regression. The result obtaining the highest R 2 was chosen as the adapted curve for the expansion trajectory. Thus, there were three patterns for urban expansion: (1) exponential expansion (14 cities); (2) expansion (10 cities); and (3) s-shaped expansion (12 cities). Exponential expansion indicates that the urban growth rate increased significantly in recent years; linear expansion indicates a relatively stable development pace; and s-shaped expansion can be explained as a city experiencing an initial slow growth followed by a rapid growth finally turning into a stage with slower development. It should be noted that, in the literature, the urbanization curve is generally based on the urban population, and the s-shaped curve has attracted the most attention. Although the urban expansion curve in this research was generated by urban area, it is able to present similar information.   As for the expansion area (Figure 7), Shanghai tops the list, with 1716.03 km 2 taken up by urban expansion. For one thing, the original urban area in Shanghai was large. For another thing, the Pudong district, which developed from the late 20th century and is regarded as the second national new district after Shenzhen, promoted the urban expansion in Shanghai. Meanwhile, Lhasa, restricted by natural conditions, shows the smallest increase in urban extent, with only 96.42 km 2 taken up by urban expansion. Overall, China has experienced dramatic urban expansion; however, influenced by a development basis, government policies, and natural conditions, the urban expansion shows great heterogeneity across cities.   To further explore the expansion patterns for different cities, a map illustrating the urban expansion trajectories is provided in Figure 8. For the city areas generated in the different study years, a linear function (linear expansion), exponential function (exponential expansion), and logistic function (s-shaped expansion) were applied to fit the regression. The result obtaining the highest R 2 was chosen as the adapted curve for the expansion trajectory. Thus, there were three patterns for urban expansion: (1) exponential expansion (14 cities); (2) expansion (10 cities); and (3) s-shaped expansion (12 cities). Exponential expansion indicates that the urban growth rate increased significantly in recent years; linear expansion indicates a relatively stable development pace; and s-shaped expansion can be explained as a city experiencing an initial slow growth followed by a rapid growth finally turning into a stage with slower development. It should be noted that, in the literature, the urbanization curve is generally based on the urban population, and the s-shaped curve has attracted the most attention. Although the urban expansion curve in this research was generated by urban area, it is able to present similar information.
Interestingly, there are some connections between the growth patterns and the cities' comprehensive competitiveness. Referring to the report in 2016 by China Business News Weekly, the 36 major cities were classified into five levels: first-tier cities, new first-tier cities, second-tier cities, third-tier cities, and fourth-tier cities (Figure 7). From Figure 8, it can be seen that all the third-tier and fourth-tier cities correspond to exponential expansion, suggesting that the development level of these cities is lower as they are still in the rapid expansion stage. All the first-tier cities, except for Beijing, have experienced s-shaped expansion, indicating that these cities tend to hold a stable urban layout after high-speed development, so that the urban growth rate then slows down. Differing from the cities mentioned above, the urban expansion trajectories for the new first-tier and second-tier cities present three patterns. Specifically, most of the second-tier cities show a linear expansion trend, demonstrating that the development of these cities has been relatively steady, without undergoing high-speed expansion. Half of the new first-tier cities have experienced rapid development in recent years, especially for those cities where new districts have been constructed, such as Tianjin Binhai New District and Qingdao Huangdao District. It can therefore be recognized that urban planning policies can greatly affect development progress.
first-tier and second-tier cities present three patterns. Specifically, most of the second-tier cities show a linear expansion trend, demonstrating that the development of these cities has been relatively steady, without undergoing high-speed expansion. Half of the new first-tier cities have experienced rapid development in recent years, especially for those cities where new districts have been constructed, such as Tianjin Binhai New District and Qingdao Huangdao District. It can therefore be recognized that urban planning policies can greatly affect development progress.

Impervious Surface Variation in the Old City Zone
The urban expansion rates and modes lend support to the urban expansion laws, while the corresponding internal evolution laws remain unknown. To better describe the variation inside the urban area, we defined the urban area in the first study year as the target old city zone, and we tried to explain the internal structure variation in this zone by ISA change.
We mapped the evolution of ISA in the old city zones of 36 major cities during the past three decades (Figure 9), and an interesting trend was found: 31 cities showed the tendency of ISA first increasing and then decreasing. In general, natural land is converted into built-up land, corresponding to the progress of urbanization, so it is easy to understand the initial increase of ISA. However, how can we explain the subsequent decrease in ISA? This phenomenon may be rooted in the widespread demolition of old buildings. In China's cities, a common phenomenon is to dismantle narrow built-up areas to build new housing estates or parks with high greening rates, resulting in a decrease of the ISA. Another major factor accounting for this phenomenon may be the

Impervious Surface Variation in the Old City Zone
The urban expansion rates and modes lend support to the urban expansion laws, while the corresponding internal evolution laws remain unknown. To better describe the variation inside the urban area, we defined the urban area in the first study year as the target old city zone, and we tried to explain the internal structure variation in this zone by ISA change.
We mapped the evolution of ISA in the old city zones of 36 major cities during the past three decades (Figure 9), and an interesting trend was found: 31 cities showed the tendency of ISA first increasing and then decreasing. In general, natural land is converted into built-up land, corresponding to the progress of urbanization, so it is easy to understand the initial increase of ISA. However, how can we explain the subsequent decrease in ISA? This phenomenon may be rooted in the widespread demolition of old buildings. In China's cities, a common phenomenon is to dismantle narrow built-up areas to build new housing estates or parks with high greening rates, resulting in a decrease of the ISA. Another major factor accounting for this phenomenon may be the reconstruction of urban villages. Urban villages are considered as a product of the traditional urban-rural dual system. Specifically, the rural land belongs to "collective ownership", while urban land belongs to "national ownership". Thus, the profit-oriented land expropriation policies initiated during city development tend to absorb farmland instead of rural residential areas, resulting in rural areas becoming surrounded by an urban area, which is a phenomenon that is known as the "urban village". These areas are likely to be unregulated, so that the land distribution in these areas is commonly unreasonable with compact houses, narrow roads, and limited green lands (Figure 10a,c), greatly affecting the internal layout in an urban area. Therefore, the reconstruction of urban villages is aimed at adjusting the land structure, especially dismantling and reconstructing compact buildings. In addition, affected by the poor quality of the construction and unreasonable real estate development, buildings in an urban area generally show a short life cycle. It has been reported that the average building life cycle in China is only 30 years [77]. Large parts of demolished areas are transformed into modern communities (Figure 8a,b) and the other parts are transformed into public places (Figure 10c,d), because of the increasing need for an improved living environment. All these factors lead to the decrease of ISA in the old city zone.
buildings. In addition, affected by the poor quality of the construction and unreasonable real estate development, buildings in an urban area generally show a short life cycle. It has been reported that the average building life cycle in China is only 30 years [77]. Large parts of demolished areas are transformed into modern communities (Figure 8a,b) and the other parts are transformed into public places (Figure 10c,d), because of the increasing need for an improved living environment. All these factors lead to the decrease of ISA in the old city zone.  Except for the trend of first increasing and then decreasing, there are also some other evolution trends for the ISA in a few cities. The ISA in the old city zones of Ningbo and Changchun present a trend of steady decrease. The reason why these cities did not show an initial stage of increasing ISA may be that the built-up areas in the old city zone of the two cities were relatively saturated by 1986. The trend in Shijiazhuang and Yinchuan was an initial slight decrease followed by an increase. Compared with the aforementioned cities, the increasing trend in recent years has been significant, which is probably attributable to government decisions. For Shijiazhuang, the increasing trend may be a result of the "great changes within three years" policy that has been enacted since 2007. For Yinchuan, the reason for the increasing trend may be the series of giant building projects undertaken Except for the trend of first increasing and then decreasing, there are also some other evolution trends for the ISA in a few cities. The ISA in the old city zones of Ningbo and Changchun present a trend of steady decrease. The reason why these cities did not show an initial stage of increasing ISA may be that the built-up areas in the old city zone of the two cities were relatively saturated by 1986. The trend in Shijiazhuang and Yinchuan was an initial slight decrease followed by an increase. Compared with the aforementioned cities, the increasing trend in recent years has been significant, which is probably attributable to government decisions. For Shijiazhuang, the increasing trend may be a result of the "great changes within three years" policy that has been enacted since 2007. For Yinchuan, the reason for the increasing trend may be the series of giant building projects undertaken since 2008 (e.g., Ningxia Software Park). The increasing trend of ISA in the two cities indicates that the city construction in recent years has taken place in undeveloped areas or green land. In particular, Haikou is the only city that has shown a tendency for consistent increase. As a capital city that was set up in 1988, Haikou's foundation has been limited by the municipal facilities. As a result, a great deal of bare land in the old city zone has been converted into ISA during the urbanization progress, and the ISA has not yet reached a state of complete saturation in the old city zone.

Urban Expansion Interaction with Economy and Population
It has been reported that urban expansion is highly correlated with the economy and population growth [44,[78][79][80][81][82]. In this research, we conducted a regression analysis using Gross Domestic Product (GDP) per capita (GDPPC) to represent the economy, and total population at the year-end (TP) for the population. The Ordinary Least Square (OLS) regressions between ISA and GDPPC, and between ISA and TP showed that the determination coefficients (R 2 ) for most cities were greater than 0.9 (Table S2). The smallest determination coefficients were greater than 0.84 (for GDPPC) and 0.69 (for TP), verifying the high correlation between urban expansion and economy, and between urban expansion and population in China's cities. To further reveal the interaction, the potential driving mechanism was investigated with the help of GCT [59,83,84] (Supplementary Information Text S3). Since the ISA used in this paper were not time-continuous, a regression between the year and the ISA was conducted to obtain the probable ISA in missing years. In detail, the ISA variation was fitted by linear, exponential and s-shaped functions, and the best fitting function was selected. Thus, the interpolation could be carried out following the selected function and a time-continuous series of ISA could be gained. Next, the Augmented Dickey-Fuller (ADF) test was applied to the time-continuous series of ISA, GDPPC, and TP. Co-integration was found in all cities except for Kunming by the Johansen co-integration test. Therefore, the GCT was applied to 35 (Kunming excluded) of the 36 cities. Furthermore, since the Granger causality result may not be the same among different lag orders, the result with the optimal lag order (verified by the Akaike information criterion and the Schwarz information criterion) was chosen. If there was no granger causality when using the optimal lag order, the result found in most of the lag orders was chosen. Based on the GCT, if event a "Granger-causes" event B, event a is taken as a potential driver.
From the GCT in the 35 cities (Table S3; Figure 11), there existed an obvious north-south differentiation for the interaction between ISA and TP. Taking the Yangtze River as the division between north and south China, most of the cities in the south indicated that TP Granger-causes ISA. Meanwhile, in the north, most of the cities indicated that ISA Granger-causes TP. a possible explanation for this may be that, for the south, it is a rapidly developed zone that has benefited from earlier reform and opening-up. Therefore, a number of young adults are continually being attracted to seek work opportunities in the south, bringing new vigor and vitality into the urban construction. As a result, the population has become the major driver of the southern cities' development. An exception is the Zhejiang province, where many people are willing to travel for business, likely resulting in ISA Granger-causing TP for Hangzhou and Ningbo. For the north, in contrast, the opening-up took place later. Thus, many people in the north have rushed to the south for work opportunities. Hence, most of the northern cities should conduct urban construction first, then attract the backflow of population. ISA has therefore become the major driver of the northern cities' development. Remote Sens. 2019, 11, x FOR PEER REVIEW 13 of 19 Figure 11. Granger causality types in China's 36 major cities. The blue line represents the Yangtze River.

The Differences in Urban Expansion Trajectories
In our study, we found three different growth patterns for urban expansion in China, i.e., s-shaped, linear, and exponential. According to the current research, some researchers insist that the s-shaped curve is more suitable for developed countries [88,89], while a j-shaped curve (similar to the exponential curve) is more appropriate for developing countries [90]. However, our results have indicated that great heterogeneity exists for urban expansion trajectories, even in the same country.
Moreover, although the urban expansion trajectories are relatively uniformed for cities in the same development level, there are still some exceptions. For example, as a first-tier city, Beijing demonstrates a linear expansion trend rather than a s-shaped trend because it developed in much earlier years and might have reached a stable development stage in the 1980s. Haikou and Urumqi, as the cities list at the end of second-tier cities, present urban expansion trajectories more similar with third-tier cities, i.e., exponential growth. We can therefore recognize that those cities also need rapid urban expansion to develop themselves. In addition, Chengdu as the top city of new first-tier cities, shows a similar s-shaped trend of urban expansion trajectory with first-tier cities. We can also infer that the urban growth in Chengdu has reached a stable stage, and it may get closer to the development of first-tier cities in the future.
In general, a fact emerged that both linear and exponential trends represented initial stages of s-shaped growth. In other words, for those cities without a s-shaped trend, we suggested that they would turn into a s-shaped development in the future. In reverse, for a city in China, when its growth pattern was figured out (exponential, linear or s-shaped), we could roughly point out its current development level.

The Limitations of Our Work
Even though we have tried to take all aspects of theentire research matter into consideration, there are some limitations to our work. Firstly, on account of the effects of cloud cover and sensor For the interaction between ISA and GDPPC (Table S3; Figure 11), the majority of the cities (27 cities) indicated that ISA Granger-causes GDPPC, confirming that the economic growth in China's cities has been greatly dependent on ISA expansion in the past three decades [85][86][87]. It is frustrating that some cities in China have expanded blindly to pursue economic growth. In particular, the GDP growth is heavily reliant on real estate in some cities, resulting in wild urban space expansion, bringing some negative effects. To solve these problems, the government has promulgated policies for controlling urban size since 2013 ("National Plan on New Urbanization 2014-2020"). As expected, only a few cities (eight cities) indicated that GDPPC Granger-causes ISA. All these cities are located to the north of the Yangtze River. Supported by resource superiority (Dalian, Taiyuan, Urumqi) or policy investment (Shijiazhuang, Lanzhou, Xining, Wuhan, Hefei), the economic returns of these cities are mainly used for urban construction.
However, it should be recognized that GCT does not prove that one variable "causes" another except in the sense that it precedes the other in time. Even though one factor is likely to promote the other when one factor happened earlier than the other, in the future, more concerns about land requisitions and government policies should be involved to assist the GCT results in pointingto a more "real" causality.

The Differences in Urban Expansion Trajectories
In our study, we found three different growth patterns for urban expansion in China, i.e., s-shaped, linear, and exponential. According to the current research, some researchers insist that the s-shaped curve is more suitable for developed countries [88,89], while a j-shaped curve (similar to the exponential curve) is more appropriate for developing countries [90]. However, our results have indicated that great heterogeneity exists for urban expansion trajectories, even in the same country.
Moreover, although the urban expansion trajectories are relatively uniformed for cities in the same development level, there are still some exceptions. For example, as a first-tier city, Beijing demonstrates a linear expansion trend rather than a s-shaped trend because it developed in much earlier years and might have reached a stable development stage in the 1980s. Haikou and Urumqi, as the cities list at the end of second-tier cities, present urban expansion trajectories more similar with third-tier cities, i.e., exponential growth. We can therefore recognize that those cities also need rapid urban expansion to develop themselves. In addition, Chengdu as the top city of new first-tier cities, shows a similar s-shaped trend of urban expansion trajectory with first-tier cities. We can also infer that the urban growth in Chengdu has reached a stable stage, and it may get closer to the development of first-tier cities in the future.
In general, a fact emerged that both linear and exponential trends represented initial stages of s-shaped growth. In other words, for those cities without a s-shaped trend, we suggested that they would turn into a s-shaped development in the future. In reverse, for a city in China, when its growth pattern was figured out (exponential, linear or s-shaped), we could roughly point out its current development level.

The Limitations of Our Work
Even though we have tried to take all aspects of the entire research matter into consideration, there are some limitations to our work. Firstly, on account of the effects of cloud cover and sensor malfunction, we could not obtain high-quality Landsat remote sensing images of 36 cities for all the years from 1986 to 2015. To solve this problem, we tried to recover the missing regions in the Landsat-7 images using a gap-filling method, and used temporal interpolation in the GCT test, which may have brought about some uncertainties which were difficult to quantify. Secondly, the population and GDP statistical data used in this paper were counted at the city level, but the ISA was computed on the main urban areas. Nevertheless, we considered this slight mismatch to be acceptable because the population and economic earnings were mostly concentrated in these areas. However, if district-level statistical data were available, the results would be more reliable. Thirdly, although the sub-pixel unmixing technique was used to extract the impervious surfaces, some bias inevitably exists. This is also one of the reasons for the small fluctuations in the impervious surface evolution for the old city zones shown in Figure 9. If the errors could be further reduced, we could make a more refined analysis of the urbanization progress by considering more influencing factors.

Conclusions
In this paper, we mapped and analyzed three decades of urban expansion in China's 36 major cities. Both homogeneity and heterogeneity were established. With regard to homogeneity, almost all of China's major cities have experienced rapid urban expansion. Meanwhile, a remarkable trend was found in the ISA variation in the city cores of the majority of China's cities, i.e., ISA first increasing and then decreasing. This can be attributed to the urban development progress in China, especially the reconstruction of urban villages. Furthermore, the interaction between ISA and the economy tends to be ISA promoting the economy, confirming that the Chinese economy is highly dependent on urban expansion. In terms of heterogeneity, constrained by development basis, government policies, and natural conditions, the urban expansion trajectories show three modes: exponential, linear and s-shaped. Interestingly, there exists an obvious north-south territorial differentiation for the interaction between ISA and population. Cities in the north tend to indicate that ISA Granger-causes population, while cities in the south are more likely to indicate that population Granger-causes ISA. This phenomenon is in accordance with the fact that a lot of young adults from the north seek work opportunities in the south.
The results suggested that the policies issued by both central and local government have had a great impact on the progress of urbanization in China [58,85]. For instance, since Shenzhen was set up as a special economic zone in 1980, it has changed from a small fishing village into an international metropolis, with a GDP that transcends that of Hong Kong. Furthermore, since the Pudong New District was developed in 1990, economic efficiency has increased 118 times in 2014 compared to the initial period. In 2017, the central government decided to develop Xiong'an New District and tried to build it into a national new district after Shenzhen and Pudong, aiming to create a new era for reform and opening-up. The positive effects exerted by the national planning of new districts have been practically validated. However, excessively pursuing local GDP growth in some local policies has resulted in some blind expansion for new district development, which brings about not only traditional "city diseases", but also some new negative effects. For example, some uninhabited new district, the so-called "ghost towns", have resulted in a great waste of resources. In addition, the fact that economic earnings and government receipts depend more on real estate may lead to weak growth in the subsequent development. Even so, "creating another new town" is still the development goal of many cities. To deal with these problems, the central government has introduced policies (e.g., "New Urbanization Planning" in 2014 [91]) to strictly control urban size. Accordingly, local government should abandon the blind pursuit of political achievement and actively adjust city development strategies, insisting on sustainable development.