Automatically Monitoring Impervious Surfaces Using Spectral Generalization and Time Series Landsat Imagery from 1985 to 2020 in the Yangtze River Delta

Accurately monitoring the spatiotemporal dynamics of impervious surfaces is very important for understanding the process of urbanization. However, the complicated makeup and spectral heterogeneity of impervious surfaces create difficulties for impervious surface monitoring. In this study, we propose an automatic method to capture the spatiotemporal expansion of impervious surfaces using spectral generalization and time series Landsat imagery. First, the multitemporal compositing and relative radiometric normalization methods were used to extract phenological information and ensure spectral consistency between reference imagery and monitored imagery. Second, we automatically derived training samples from the prior MSMT_ IS30-2020 impervious surface products and migrated the surface reflectance of impervious surfaces in the reference period of 2020 to other periods (1985–2015). Third, the random forest classification method, trained using the migrated surface reflectance of impervious surfaces and pervious surface training samples at each period, was employed to extract temporally independent impervious surfaces. Further, a temporal consistency check method was applied to ensure the consistency and reliability of the monitoring results. According to qualitative and quantitative validation results, the method achieved an overall accuracy of 90.9% and kappa coefficient of 0.859 in identifying the spatiotemporal expansion of impervious surfaces and performed better in capturing the impervious surface dynamics when compared with other impervious surface datasets. Lastly, our results indicate that a rapid increase of impervious surfaces was observed in the Yangtze River Delta, and the area of impervious surfaces in 2000 and 2020 was 1.86 times and 4.76 times that of 1985, respectively. Therefore, it could be concluded that the proposed method offered a novel perspective for providing timely and accurate impervious surface dynamics.


Introduction
Impervious surfaces, mainly composed of anthropogenic materials that prevent water from penetrating into the soil [1,2], are important indicators for understanding the process of urbanization and assessing the environmental quality [3]. According to the prospects of the United Nations in 2018, 55% of the world's population lives in urban areas, and the proportion of the global urban population will reach 68% in 2050 [4]. Thus, a series of social, environmental, and climatic problems are emerging with the continuous expansion of impervious surfaces, including traffic congestion, urban heat islands, food shortages, water quality deterioration, and loss of biodiversity [5][6][7], which directly affect human well-being and health [8]. Therefore, accurately monitoring the spatiotemporal dynamics of impervious surfaces is greatly important for understanding urbanization's impacts on environmental quality and achieving the goal of sustainable urban development.
Over the past decades, impervious surface monitoring has been experiencing a transition from a coarse resolution of 1 km to a fine spatial resolution of 30 m with the development of remote sensing techniques and increasing computation and storage abilities [8][9][10][11]. Correspondingly, a series of impervious surface monitoring methods have also been proposed for various satellite data with different spatial and temporal resolutions. Overall, these monitoring methods could be grouped into two main categories, including time series monitoring methods [12,13] and independent classification methods [5,[14][15][16]. Specifically, the independent classification strategy is to split the time series into multiply independent time points and then produce these temporally independent classifications, and it can be further divided into preclassification and postclassification [13]. The preclassification method treats land cover changes as independent land cover types and then collects stable and changed training samples to directly classify multidate imagery for generating temporally consistent impervious surface results [5,17], while the postclassification uses multiepoch training samples to independently train the classifier on the corresponding imagery and then combines the multiepoch impervious surface classification results to monitor the spatiotemporal expansion of impervious surfaces [15,18,19]. However, it is time-consuming and labor-intensive when collecting the training samples, especially for preclassification methods, which directly give the strategy great difficulty in monitoring impervious surface expansions over large areas. Recently, to achieve automatically monitoring large-area impervious surface dynamics, the combination of impervious surface spectral index and thresholdsegmentation was proposed to generate multiply temporally independent classifications. Specifically, it builds a series of impervious-sensitive indexes and further uses optimal thresholds to identify the impervious surface extents period by period, which can achieve fast and automatic monitoring of impervious surface dynamics, but there are great difficulties in finding the optimal thresholds to accurately separate the impervious surfaces from pervious surfaces, especially in arid regions [1,20]. For example, Liu et al. [11] combined the proposed NUACI (Normalized Urban Areas Composite Index) and the local adaptive thresholds to produce multitemporal global impervious surface maps for the 1985-2015 period with a five-year interval, but the kappa coefficients of the products are 0.43-0.50 at the global level.
The time series monitoring strategy, taking advantage of free access to Landsat and Sentinel-2 imagery, has been widely used in monitoring the dynamics of the terrestrial environment including deforestation and afforestation [21][22][23] as well as impervious surface expansion [13,24]. The strategy has a stronger ability to make use of the continuous Landsat imagery compared with the previous strategy. However, the time series change detection monitoring method still needs independent classification results as the prior knowledge when monitoring the impervious surface expansion. For example, Liu et al. [13] excluded areas that were temporally persistent or irrelevant land cover changes by using the independent land cover classifications in the start and end years, which also required manual involvement. Then, interannual change detection methods have also been used to capture dynamic impervious surfaces [13], but there has been uncertainty in the methods, which transforms the final monitoring results. For example, the new Continuous Monitoring of Land Disturbance (COLD) dataset was validated to achieve an omission error of 27% and a commission error of 28% for monitoring land disturbance [22]. Therefore, the time series monitoring strategy also has difficulty automatically in monitoring large-area impervious surface dynamics.
The expansion of impervious surfaces is the most irreversible process among all land cover changes, especially in developing areas [13,19]. Therefore, if the maximum impervious surface extent can be determined in advance, stable pervious surface samples could be easily derived, so monitoring impervious surface dynamics could focus only on how to automatically derive multitemporal impervious surface samples. Meanwhile, to avoid the difficulty of manually collecting samples, spectral generalization has been demonstrated to perform well for automatic land cover mapping and monitoring [25][26][27]. For example, Zhang et al. [25] used spectral generalization to generate multitemporal land cover products in China and achieved an overall accuracy of 80.7%. However, Hu and Liu [28] found that the performance of spectral generalization had a positive relationship with the spectral consistency between reference imagery and unclassified imagery. Therefore, spectral generalization could be used to automatically monitor impervious surfaces after ensuring consistency between reference imagery and unclassified imagery.
The present study proposes a novel method for automatically and accurately monitoring the spatiotemporal expansion of impervious surfaces using spectral generalization and time series Landsat imagery. To achieve this goal, we migrated the surface reflectance of impervious surfaces in the reference period of 2020 to other periods and combined independent random forest models and a temporal consistency check method to generate spatiotemporal expansion results of impervious surfaces. Lastly, the proposed method was tested in the Yangtze River Delta and assessed by validation samples and cross-comparisons. The results indicate that the proposed method accurately captured the spatiotemporal expansion of impervious surfaces automatically.

Study Areas and Datasets
2.1. Study Areas and Time Series Landsat Imagery. The Yangtze River Delta (YRD) is one of the fastest economic growth areas since China opened its doors to the world in the late 1970s [5], with some of the most active development of China's economy and the highest level of openness and innovation. Located at the Yangtze River and Donghai convergence, it is an important triangle-shaped megalopolis composed of Shanghai, Zhejiang, Jiangsu, and Anhui provinces, with an area of 35,800 km 2 ( Figure 1). The study area has a population of approximately 227 million and a total gross domestic product of 23.72 trillion yuan at the end of 2019 accounting for nearly one-fourth of the country (National Development and Reform Commission). The fast economic development and population growth have led to a sharp urbanization process during the past decades causing drastic land cover changes, mainly including the change from agricultural and natural land cover types to impervious surfaces.

2
Journal of Remote Sensing To completely monitor the spatiotemporal expansion of impervious surfaces, all available Landsat 5, 7, and 8 imagery from 1984 to 2020 was collected via Google Earth Engine. For each image, six optical reflectance bands were retained along with one quality assessment band, which labeled "bad quality" observations including snow, cloud, cloud shadow, and saturated pixels using the CFmask algorithm [29]. In order to minimize atmospheric effects, all Landsat 8 images were atmospherically corrected using the Land Surface Reflectance Code (LaSRC) method [30], and the Landsat 5 and 7 images were corrected using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm [31] by the United States Geological Survey and then archived on the Google Earth Engine platform. In total, 27,973 Landsat images were used, including 4,880 Landsat 8 OLI, 11,074 Landsat 7 ETM+, and 12,019 Landsat 5 TM images. Figure 2 illustrates the spatial distribution of availability of time series Landsat imagery after masking the "bad quality" observations over the period of 1985-2020 with the interval of 5 years. Intuitively, more Landsat imagery is available recently, because later Landsat satellites had higher capacities for onboard recording.

Prior Impervious Surface Products in 2020.
To automatically monitor the spatiotemporal expansion of impervious surfaces, prior impervious surface data from 2020 was necessary. In this study, MSMT_IS30-2020 (global 30 m impervious surface product using multisource and multitemporal datasets in 2020) was used as the reference dataset for deriving the training samples and to determine the maximum impervious extent when monitoring the impervious surface dynamics in Section 3. The dataset was produced by combining time series Sentinel-1 SAR, Landsat 8 surface reflectance, and SRTM topographical imagery and validated to achieve higher accuracy than other global 30 m impervious surface products with an overall accuracy of 95.1% and kappa coefficient of 0.898 [1].

Validation Dataset.
To quantitatively analyze the performance of time series impervious surface products, a total of 3,845 validation samples were collected, including 2,160 pervious surface samples and 1,685 impervious surface samples covering a long-term time series from 1985 to 2020 ( Figure 3). Specifically, as Google Earth provided time series high-resolution remote sensing imagery, each validation sample was labeled as "pervious surface" or "specific change year" for the impervious surface. It should be noted that the time series Landsat images from 1984 to 2000 were also imported to assist the visual interpretation because highresolution images on Google Earth were relatively sparse before 2000. In addition, it was fair to validate our time series products with 30 m spatial size in the high-resolution imagery [20], because the spatial structures of impervious surfaces were broken compared with the natural surfaces. If the impervious area in the 30 m × 30 m window was less than the threshold of 50%, the validation sample was directly discarded from the validation dataset. Furthermore, to minimize the effect of geometry registration between validation samples and our products and further improve the confidence of validation samples, the geolocations of these rural impervious surface samples which located in the boundary  Figure 1: The geographical location of the Yangtze River Delta. The background image was the Landsat maximum NDVI composited imagery in 2020.
3 Journal of Remote Sensing of the impervious objects (such as buildings and roads) were positioned in the center of the objects. It should be noted that a total of 37 rural impervious surface samples have been moved after our statistics. The multitemporal impervious surface validation dataset in the Yangtze River Delta is free to access at doi:10.5281/zenodo.4905198.

Method
In this study, we assumed that the process of transforming pervious surfaces into impervious surfaces was irreversible during the period 1985 to 2020. Therefore, the pervious surface samples derived from the prior impervious surface products in 2020 would be kept stable during the period. As for the impervious surface samples, collecting the corresponding samples at each period would require manual participation, so we proposed to temporally generalize the reflectance spectra of impervious surfaces in 2020 to other periods of 1985 to 2015 to automatically monitor the spatiotemporal expansion of impervious surfaces. Intuitively, Figure 4 illustrates the flowchart of the proposed method, which mainly includes two stages: deriving training reflectance spectra and spectral generalization for monitoring time series impervious surfaces. Detailed procedures of the proposed method are described below.
3.1. Deriving Training Reflectance Spectra. It was impractical to collect timely and consistent training samples to monitor spatiotemporal expansion of impervious surfaces, especially for the immense study area. Therefore, the spectral generalization method was proposed to automatically migrate the reflectance spectra of impervious surfaces in the reference period (2020) to other periods . It should be noted that only the reflectance spectra of impervious surfaces needed to be migrated because the pervious surface samples derived from the MSMT_IS30-2020 products were kept stable during the whole period of 1985-2020 based on the irreversibility assumption. Namely, the reflectance spectra of pervious surfaces at other periods were extracted from the seasonal composites acquired at the corresponding period using the same pervious samples at 2020. The keys of this procedure were to ensure the spectral consistency between the reference period and target periods and deriving training samples based on the MSMT_IS30-2020 products.

Time Series Landsat Preprocessing. The time series
Landsat images used in this study have been atmospherically corrected to the surface reflectance, but the cloud, shadow, and snow pixels were invalid observations for monitoring the impervious surfaces, so all poor-quality observations in the Landsat imagery were first masked using the corresponding quality assessment band.

Journal of Remote Sensing
Many studies have demonstrated that multitemporal information has a positive contribution when separating impervious surfaces from other pervious surfaces (e.g., cropland, water, bare land, and forest) [1,32,33]. For example, our previous study [1] found that there was a large degree of complementarity between images from different seasons. To ensure the accuracy of the final time series impervious surface products in this study, it was necessary to composite the time series Landsat imagery into multitemporal variables. According to the reviews of Gómez et al. [34], there were mainly two options to capture multitemporal information from time series Landsat imagery: pixel-based image compositing and disentangling spectrotemporal feature space into different variables. It should be noted that the former selects the best available observation subject from Landsat imagery based on user-defined criteria so the composited imagery still contains the characteristics of surface reflectance [35]. Conversely, the latter option transforms the spectrotemporal feature space into other variables. For example, the multitemporal spectral data can be summarized into statistical metrics (e.g., average, minimum, and maximum) [36] and change metrics (magnitude, duration, and slope) [34], so the transformed variables cannot be equal to the original Landsat surface reflectance images.
In this study, the spectral generalization was to migrate the reflectance spectra of impervious surfaces in the reference period to other periods; therefore, the pixel-based image compositing strategy was selected to capture multitemporal information from time series Landsat imagery. Specifically, the best-available-pixel (BAP) compositing approach was used to generate cloud-free surface reflectance composites, calculating four scores (sensor score, day of year score, distance to cloud score or cloud shadow score, and opacity score) for each pixel and finding the best observation with the maximum comprehensive total score [37]. Furthermore, to capture intra-annual variability from the time series Landsat imagery, the BAP method was applied for each season so the multiple seasonal surface reflectance imagery (spring, summer, autumn, and winter) was generated over each period. The six spectral reflectance bands were imported to the multitemporal feature space for each BAP composited image along with the normalized difference vegetation index (NDVI), normalized difference water index (NDWI), and normalized difference built-up index (NDBI), which have been demonstrated to contribute to separating impervious surfaces from natural pervious surfaces [1,11]. This provided a total of 36 features for four seasonal surface reflectance composites.
Our previous study demonstrated that the spectral consistency between reference imagery and target imagery had an influence on the accuracy of the spectral generalization method [28,38]. In this study, the BAP seasonal compositing method tried to guarantee the phenological consistency of seasonal composites from multiple periods, but due to differences in the spectral response function (OLI, ETM+, and TM) and atmospheric correction methods (LEDAPS and LaSRC) as well as the temporal differences between reference imagery and classified imagery, the relative radiometric normalization method should be applied to the BAP seasonal composites before deriving training samples. Specifically, as the spectral generalization migrated the reflectance spectra of impervious surfaces in 2020 to other periods, the seasonal reflectance composites in 2020 (ρ 2020 ðλ i , x, yÞ) and other periods (ρ T ðλ i , x, yÞ) were the dependent variables and independent variables, respectively, where α T,i and β T,i are the slope and intercept of the linear regression formula for relative normalization between the reflectance band i at the period of T with the reference period of 2020.

Deriving Training Samples from the MSMT_IS30-2020
Products. Many studies have explained that the representativeness and reliability of training samples have great influence on classification accuracy [39,40]. In this study, MSMT_IS30-2020 was validated to achieve an overall accuracy of 95.1% and kappa coefficient of 0.898, so the training samples derived from the MSMT_IS30-2020 products had enough confidence for monitoring impervious surface expansion. The pervious samples derived from 2020 could be directly used in the other periods based on the irreversibility assumption. As for the representativeness of training samples, the distribution, size, and balance of training samples were all important factors. For example, Jin et al. [41] found that training samples with equal distribution (the number of training samples was the same for all land cover types) achieved higher producer's accuracy than that of proportional distribution for identifying urban and nonurban land cover types. In this study, the focus was to monitor the spatiotemporal expansion of impervious surfaces, a sparser land cover type compared with the pervious surfaces. Therefore, training samples with equal distribution were used to ensure the rationality of training samples and capture all relevant spectral heterogeneity within impervious surfaces. For the training sample size, Zhu et al. [42] quantitatively analyzed the relationship between sample size and classification accuracy and suggested a maximum of 8,000 training samples for the dominant land cover type. Therefore, we randomly generated the 16,000 training samples based on the MSMT_IS30-2020 products.

Spectral Generalization for Monitoring Spatiotemporal
Expansions of Impervious Surfaces. This study differed from traditional spectral generalization methods to migrate the reflectance spectra of all land cover types in that we only migrated the reflectance spectra of impervious surfaces in 2020 to other periods according to the irreversibility assumption. Afterwards, the classification models were separately trained by combining the migrated reflectance spectra of impervious surfaces with the training spectra of pervious surfaces at different times. We trained seven independent classification models to classify the corresponding impervious surface extent from 1985 to 2015.
As for the selection of a specific classifier, the random forest classification model had obvious advantages over the 6 Journal of Remote Sensing classification accuracy and computation efficiency compared with other widely used classifiers, including support vector machine, artificial neural network, and decision tree, based on our previous investigations [25,43]. It had significant robustness for handling high-dimensional data [44] and was insensitive to the low training error and training feature selection [45]. Therefore, it was used to produce the time series impervious surface products. In addition, the random forest classifier had only two adjusted parameters including the number of decision trees (Ntree) and the number of features for splitting each tree node (Mtry). The best choice for Mtry was the default value (square root of the total input features) according to comprehensive analyses in previous studies [44,46], and 500 was suggested for the Ntree because the classification error stabilized before reaching the threshold. Therefore, the default values of these two parameters were implemented to train the random forest classifier. It was necessary to use a temporal postprocessing method to ensure the consistency and reliability of the monitoring results because our time series impervious surface products were generated by independent classifications, and we assumed that the process of transforming pervious surfaces into impervious surfaces was irreversible. The "temporal consistency check" method proposed by Li et al. [18], consisting of temporal filtering and logical reasoning, was applied because it simultaneously took into account the temporal adjacent contexts and prior knowledge to suppress the illogical conversions caused by the classification errors. The "temporal filtering" procedure was used to iteratively remove the misclassification caused in the individual classifications, and "logical reasoning" was applied to ensure the logical conversion from pervious surfaces to impervious surfaces [18].
3.3. Accuracy Assessment. The "pixel-based" and "crosscomparison-based" methods were used to comprehensively assess the performance of the time series impervious surface products. The "pixel-based" method used the multitemporal validation samples to develop the confusion matrix and calculate four typical accuracy metrics, including producer's accuracy, user's accuracy, overall accuracy, and kappa coefficient [47]. It should be noted that the time series impervious surfaces products and the validation samples both contained temporal information, so we directly evaluated the accuracies of change detection performance, which treated the expansion periods as evaluation objects.
The "cross-comparison-based" method was implemented to intuitively compare our time series impervious surface results with the existing multitemporal products. Based on the review of [48], there were eight global 30 m impervious surface products for free access, but only 5 of them contained multitemporal characteristics, including the Global Human Settlement Layer (GHSL, four epochs of 1975,1990,2000, and 2014) [49], the GlobeLand30 impervious layer (three epochs of 2000, 2010, and 2020) [50], Global Artificial Impervious Area (GAIA, 34 epochs from 1985 to 2018) [8], the global multitemporal urban land (NUACI, seven epochs of 1980,1990,1995,2000,2005,2010, and 2015) [11], and the Global Urban Dynamic (GUD, 31 epochs from 1985 to 2015) [51]. It should be noted that the GAIA and GUD annual impervious surface products were grouped into eight epochs at an interval of five years for the sake of cross-comparisons.

Spatiotemporal
Expansion of Impervious Surfaces from 1985 to 2020. Figure 5 illustrates spatial patterns of time series impervious surfaces from 1985 to 2020 at an interval of five years. Intuitively, the Yangtze River Delta experienced rapid urban expansion in the past 35 years. In 1985 ( Figure 5(a)), the impervious surfaces were mainly concentrated in the large city centers and the links between cities were not close (the impervious surface pixels between adjacent cities were relatively sparse). After 2000 (Figures 5(d)-5(h)), the cities were obviously connected by the new impervious surfaces, especially in the areas surrounding Shanghai. The economic and population growth, mainly driven by manufacturing, has promoted the expansion of impervious surfaces in these areas.
The impervious surface area of the eight periods was calculated ( Figure 6) to quantitatively analyze the expansion speed of impervious surfaces from 1985 to 2020. The expansion rate of impervious surfaces in the Yangtze River Delta after 2000 was significantly faster than before 2000. Quantitative statistics show that the area of impervious surfaces in 2000 and 2020 was 1.86 and 4.76 times that of 1985, respectively. The average growth rate before 2000 was about 442.9 km 2 per year, and the growth rate after 2000 was 1,151.1 km 2 per year. This phenomenon is apparent in Figure 5; the cities expanded relatively slowly from 1985 to 2000 and significantly accelerated after 2000. Similarly, Gao et al. [5] also found that many cities in the Yangtze RD obviously accelerated their development after 2000 mainly because the industrial economy became a dominant sector in this region.
The eight-period impervious surface mapping results were composited into the spatiotemporal expansion results in Figure 7 to intuitively understand the spatiotemporal expansion of impervious surfaces from 1985 to 2020 in the Yangtze River Delta. The expansion pattern followed a diffusion-coalescence model [52], where early expansion of impervious surfaces occurred in isolated patches (diffusion), and subsequent expansions eventually filled the gaps between urban centers and connected previously isolated urban areas by new impervious areas (coalescence). From the perspective of expansion intensity, cities along the Yangtze River and coastal cities have experienced the most dramatic impervious surface expansion over the past 35 years, especially around Shanghai and Hangzhou Bay. The enlargements in Shanghai and surrounding areas are illustrated as an example in Figures 7(b)-7(d). They have formed an obvious urban agglomeration from the previous "island" cities, and a large amount of croplands have been converted into impervious surfaces, which also dominated the land cover changes around Shanghai from 1985 to 2020.   Journal of Remote Sensing and producer's accuracy of 93.9% because the old cities before 1985 occupied a larger proportion than the subsequent growth and also had relatively uniform spectral reflectance characteristics (Figure 7(c)). Next, the producer's and user's accuracies of the increased impervious surfaces in the seven periods were approximately 70%, mainly because there were misclassifications between three temporally adjacent periods. The increased impervious surfaces during 1995 to 2000 were mainly confused with the previous period of 1990 to 1995 and the next period of 2000 to 2005.
To further clarify the confusion of pervious surfaces and various periods' impervious surfaces, the confusion proportions between pervious surfaces and impervious surfaces in Table 1 are calculated and illustrated in Figure 8. The pervious surfaces and impervious surfaces before 1985 had low confusion proportions with other types, and confusions were  9 Journal of Remote Sensing mainly concentrated in the increased impervious surfaces of seven periods. Approximately 15% of increased impervious surfaces during 1985 to 1990 were misclassified as those before 1985, and more than 10% of increased impervious surfaces during 2015 to 2020 were misclassified as those between 2010 and 2015. Generally, these confusions were usually reasonable because the conversion from natural surfaces to impervious surfaces in 30 m Landsat imagery was a gradual process, which also caused the validation samples to have temporal deviation.

Comparison with Other Existing
Products. In addition to the quantitative accuracy assessment, a qualitative crosscomparison between our results (MSMT_IS30) and other impervious surface products (including GAIA, NUACI, Glo-beLand30, GHSL, and GUD) in Shanghai, Suzhou, Hang-zhou, and Hefei is illustrated in Figure 9. Overall, MSMT_ IS30, GAIA, and GlobeLand30 captured the spatiotemporal expansion of impervious surfaces in these regions, while NUACI and GHSL suffered from overestimation in the earlier years (e.g., 1990) compared to Landsat imagery in the corresponding year. Similarly, Gong et al. [8] also found that GHSL was overestimated in the early epoch of 1990 over Beijing.
The comparisons in Shanghai showed that the NUACI and GHSL results failed to capture the expansions of Shanghai. For example, the impervious extents of NUACI and GHSL before 2000 (green color) were much larger than the actual distributions in 2000 and the increased impervious surfaces after 2000 only accounted for a small percentage, especially in the NUACI results. GlobeLand30 only provided three epochs (2000, 2010, and 2020) of impervious surface  Figure 8: The confusion proportions between pervious surfaces and impervious surfaces from various periods.

10
Journal of Remote Sensing dynamics, but we found that the impervious surface area of GlobeLand30 was obviously larger than the actual areas in Landsat imagery. For example, farmlands at peripheral cities were misclassified as impervious surfaces, and the fragmented structural features of Shanghai were also lost in Globe-Land30. The overestimation problems mainly came from the minimum 4 × 4 mapping unit for impervious surfaces in GlobeLand30 [50]. Although GAIA reflected the continuous expansion of Shanghai from the center to the surroundings, the growth ratio of the impervious surfaces in GAIA was too rapid during 2010 to 2020 compared to the actual expansion patterns. The rate of urban expansion between 2000 and 2010 was significantly greater than that of 2010 to 2020 according to the multitemporal Landsat imagery. Therefore, it could be concluded that GAIA had an overestimation problem during 2010 to 2020. Lastly, compared with the other five impervious results, there was the highest consistency between MSMT_IS30 and GUD dataset; both of them accurately captured the expansion pattern of impervious surfaces in Shanghai with a relatively slow expansion ratio during 1985 to 2000 and significantly accelerating after 2000. However, as the GUD mainly concentrated on monitoring the urban expansions, a lot of small impervious objects (including rural buildings) in the GUD were wrongly labeled (red rectangle).
Similar to the previous analysis in Shanghai, NUACI still could not capture the expansions of impervious surfaces and obviously overestimated them in the earlier stages (before 2000) in Suzhou, Hangzhou, and Hefei. GlobeLand30 accurately monitored the impervious surface dynamics in the three remaining regions, but it still could not reflect the spatial details of the spatiotemporal expansion of impervious surfaces and had a larger impervious surface area than the five other results. GHSL performed better in the three remaining regions than in Shanghai. The overestimation in earlier stages was improved, and the increased impervious surfaces after 2000 were also captured. However, GHSL underestimated these regions, especially in Hangzhou, where the maximum impervious surface extent of GHSL was smaller than that of the five other products. GAIA still had greater performance before 2010 and overestimated the increased impervious surfaces during 2010 to 2020 in Suzhou and Hangzhou. However, it was the opposite in Hefei because GAIA suffered from the overestimation problem in the early stage while follow-up monitoring was accurate. The GUD still had great performance for capturing the urban expansions, but a lot of rural impervious surfaces were missed especially in Hangzhou. Meanwhile, the GUD also suffered from the overestimation problem in the early stage in Hefei, namely, the impervious extent before 2000 in GUD was larger than the actual impervious extent. Therefore, based on the visual interpretation of multitemporal Landsat imagery, MSMT_IS30 still performed better than the five other products in the three remaining areas, where the increased impervious surfaces at different periods were accurately captured.

The Feasibility of Spectral Generalization for Monitoring
Impervious Surfaces. Due to the spatial and spectral heterogeneity of impervious surfaces, monitoring the spatiotemporal expansion of impervious surfaces is a labor-intensive and challenging work. In this study, the spectral generalization method, migrating the reflectance spectra of impervious surfaces in the reference period to other periods, was proposed to automatically produce time series impervious surface products. Many previous studies have demonstrated that the spectral consistency between reference imagery and classified imagery has a direct influence on the final classification accuracy for spectral generalization methods [28,53], so it was necessary to analyze the consistency between surface reflectance of impervious surfaces in the period of 2020 against corresponding surface reflectances in other periods. Specifically, we take the summer NIR composites as an example to explain the radiometric consistency ( Figure 10). It can be found that most scattered points are clustered around the 1 : 1 line, and there was significant consistency for NIR surface reflectance of impervious surfaces between the reference period of 2020 and other periods, achieving the mean determination coefficient of 0.761 and root mean square error of 0.015. Furthermore, from the perspective of regression coefficients, it could be found that the regression slope after 2000 was closer to the 1 : 1 line than that before 2000, which is mainly caused by the differences in the spectral response function (OLI, ETM+, and TM) and atmospheric correction methods (LEDAPS and LaSRC) as well as the temporal differences between reference imagery and classified imagery. It should be noted that NIR surface reflectance of impervious surface pixels between 2020 and other periods came from the BAP composited imagery instead of the normalized imagery after using Formula (1). Therefore, if the relative radiometric normalization is applied to ensure the spectral consistency between reference imagery and classified imagery, it will be feasible and rational to migrate the reflectance spectra of impervious surfaces in 2020 to the other periods.
In addition, another key element for automatically monitoring the expansion of impervious surfaces was to derive training samples from the prior MSMT_IS30-2020 products. Based on previous studies, the reliability and confidence of training samples had great influence on the final classification accuracy [42,54,55]. MSMT_IS30-2020 achieved an overall accuracy of 95.1% and kappa coefficient of 0.898 [1], so the derived training samples also achieved a high degree of confidence. As for the problem of whether the classification error of MSMT_IS30-2020 would be transmitted to the training samples, our previous study [1] quantitatively analyzed the relationship between the percentage of erroneous samples and the classification accuracy. The results indicated that the overall accuracy remained stable when the percentage of erroneous samples was within 20% [1]. Therefore, it can be concluded that the training samples derived from the MSMT_IS30 products were accurate and suitable for automatically monitoring the expansion of impervious surfaces.

Limitations and Prospects of the Proposed Method.
Although we used the proposed method to successfully monitor the spatiotemporal expansion of impervious surfaces in the Yangtze River Delta from 1985 to 2020, there are still several limitations for the automatic method. First, we assumed that the process of transforming pervious surfaces into impervious surfaces was irreversible during the period 1985 to 2020, which directly caused the method to fail to detect land cover changes from impervious surfaces to pervious surfaces and high-frequency changes in urban areas (such as reclamation and revegetation of abandoned residential land). Meanwhile, according to the review of articles on impervious surface monitoring, all global 30 m impervious surface dynamic researches (such as GAIA [8], NUACI [11,51], GHSL [56], and GlobeLand30 [50]) and 12 Journal of Remote Sensing most regional researches (such as Li et al. 2015 in Beijing [18], Zhang et al. 2016 in Pearl River Delta [19], and Gao et al. 2012 in Yangtze River Delta [5]) still use the irreversible assumption to monitor the time series impervious surface dynamics. Actually, a part of impervious surface pixels would undergo changes twice or more in a long time sequence especially in rapidly developing regions [13]; therefore, our further work should make full use of the continuous Landsat observations to assist the proposed method to monitor the high-frequency changes in urban areas.
Further, our previous study demonstrated that multisource and multitemporal features contributed to splitting impervious surfaces from natural pervious surfaces [1], but in this study, we only used the multitemporal information to monitor the spatiotemporal expansion of impervious surfaces. Although the results indicated that the multitemporal information performed well for monitoring the spatiotemporal expansion of impervious surfaces in the Yangtze River Delta, this does not mean that our method could cope with all complex scenarios, because some pervious surface pixels shared similar spectral characteristics with impervious surfaces, especially in arid and semiarid areas [1,57]. Therefore, further work should test the performance of the proposed method in arid and semiarid areas and simultaneously import time series SAR imagery, which could provide information on the structure and dielectric properties of surface materials to assist in more accurate impervious surface monitoring.
Lastly, to achieve automatically monitoring spatiotemporal dynamics of impervious surfaces, in this study, the prior MSMT_IS30 impervious surface products in 2020 were used to derive the confident training samples. However, as a universal monitoring method, the proposed method should not solidly rely on a one-source urban product. Therefore, our further work would test the performance of using other available urban products (such as GAIA, NUACI, and GHSL) for generating time series impervious surface products. Meanwhile, similar to our previous work in [1], we would further analyze the quantitative relationship between the impervious surface monitoring accuracy and the proportions of erroneous training samples.

Conclusion
In this study, a spectral generalization method was proposed to automatically monitor the spatiotemporal expansion of impervious surfaces from 1985 to 2020 in the Yangtze River Delta. We proposed to migrate the surface reflectance of impervious surfaces in 2020 to other periods and derived training samples from MSMT_IS30-2020 impervious surface products, which was the key for automatic monitoring. Then, the relative radiometric normalization method was applied to ensure consistency for surface reflectance of impervious surfaces between the reference period of 2020 and other periods. Afterwards, the random forest classifiers were trained using the migrated surface reflectance of impervious surfaces and pervious surface training samples at each period and further applied to generate time series independent impervious surface classification results. Lastly, the temporal consistency check method was applied to ensure the consistency and reliability of the monitoring results.
The spatiotemporal expansion of impervious surfaces from 1985 to 2020 was validated by 3,845 time series validation samples. The results indicated that the monitoring method achieved an overall accuracy of 90.9% and kappa coefficient of 0.859. Meanwhile, the cross-comparison with 13 Journal of Remote Sensing five existing 30 m impervious surface products indicated that our results performed better in capturing the impervious surface dynamics from 1985 to 2020. A rapid increase in impervious surfaces was observed in the Yangtze River Delta, especially since 2000. The results show that the area of impervious surfaces in 2000 and 2020 was 1.86 times and 4.76 times that of 1985, respectively. The average growth rate before 2000 was about 442.9 km 2 per year, and the growth rate after 2000 was 1,151.1 km 2 per year. Therefore, the proposed method demonstrated a strong ability to automatically capture the spatiotemporal expansion of impervious surfaces.

Data Availability
The multitemporal impervious surface validation dataset in the Yangtze River Delta is free to access at doi:10.5281/ zenodo.4905198.