A Novel Impervious Surface Extraction Method Based on Automatically Generating Training Samples From Multisource Remote Sensing Products: A Case Study of Wuhan City, China

Impervious surfaces caused by rapid urbanization affect the environment and increase the disaster risk. Currently, most articles have extracted impervious surfaces by manual participation for training samples with medium-and-high spatial resolution remote sensing images. Therefore, it is necessary to develop a new method for improving the efficiency of training sample acquisition and the accuracy of impervious surface extraction. In this article, a novel impervious surface extraction method is proposed based on automatically generating training samples from multisource remote sensing products. First, the preliminary sample area was constructed through the overlay analysis of the classification consistency area of three remote sensing products and homogeneous area detection based on Sentinel-2 images. Second, four spectral indices and digital surface model (DSM) elevation data were used for sample selection, and the pixels were further purified by variance purification calculation. Finally, by sample migration and random forest model training, impervious surfaces were extracted for other years with limited data. Wuhan city in China was selected as the study area due to a large number of interior objects for impervious surfaces. Sentinel-2 images from 2018 to 2020, three 30 m-resolution products, and DSM data in 2018 were used. The proposed method's extraction accuracies of impervious surfaces for Wuhan in 2018, 2019, and 2020 are 94.02%, 94.45%, and 93.87%, respectively. Additionally, with the resolution improved up to 10 m, the method is more conducive to distinguishing the boundary between impervious surfaces and pervious surfaces.

Abstract-Impervious surfaces caused by rapid urbanization affect the environment and increase the disaster risk. Currently, most articles have extracted impervious surfaces by manual participation for training samples with medium-and-high spatial resolution remote sensing images. Therefore, it is necessary to develop a new method for improving the efficiency of training sample acquisition and the accuracy of impervious surface extraction. In this article, a novel impervious surface extraction method is proposed based on automatically generating training samples from multisource remote sensing products. First, the preliminary sample area was constructed through the overlay analysis of the classification consistency area of three remote sensing products and homogeneous area detection based on Sentinel-2 images. Second, four spectral indices and digital surface model (DSM) elevation data were used

I. INTRODUCTION
I MPERVIOUS surfaces can prevent the infiltration of water into the soil, which is always accompanied by ground hardening, affects the green ecological environment of urban areas, aggravates urban heat islands and nonpoint source pollutants, and increases the risk, intensity, and duration of urban waterlogging with a heavy burden on drainage facilities [1], [2], [3], [4], [5], [6]. Efficient extraction and accurate estimation of impervious surfaces are of great significance for urban planning, environmental assessment, and ecological protection. Due to the rapid development of remote sensing technology, remote sensing images with an easy acquisition, wide coverage, high efficiency, and diverse data have been the main means of extracting impervious surfaces.
Impervious surface extraction is a subset of remote sensing image classification utilizing techniques, such as spectral mixture analysis [7], [8], [9], index method [10], [11], [12], regression model [13], [14], [15], decision trees [16], [17], random forest [18], [19], [20], [21], and deep learning [22], [23], [24], [25]. Among them, the extraction accuracy of impervious surfaces through spectral mixture analysis is limited by the large differences between endmember spectra and actual images. Especially, when obtaining manual endmembers from This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ image analysis, it is difficult to select good training samples for the mixed pixels classification [26], which is not conducive to automatic and rapid extraction of impervious surfaces. The index method is based on the spectral index, which is used to distinguish whether a pixel belongs to the impervious surface according to the value of the corresponding index. It adopts the form of combining multiple indices or using the index of impervious surfaces as auxiliary data for sample extraction and other classification methods. The regression model is greatly affected by the time of image acquisition. Meanwhile, the representative methods of sample selection and the minimal number of necessary samples are important to improve the efficiency of the regression model [27]. Decision trees usually need train samples and add original image bands and spectral indices as new bands to establish classification rules. The small number of training samples might lead to the loss of important information, and the classification rules could not be well explained, but a large number of training samples could lead to the increase of the nodes for the decision tree and interfere with its accuracy [28]. Random forest [29], [30] can achieve reasonable results even without hyperparameter adjustment, and the relative importance of each feature in the prediction process can be easily measured, which is simple and flexible to use. Although the random forest classifier is sensitive to the sampling design, it has been proven to be suitable for classifying hyperspectral data [31], [32]. Deep learning can integrate various features of remote sensing images, but it relies on using a complex model with poor interpretability, high hardware requirements, and many training samples. In particular, the number and diversity of training samples are essential to the multilayer neural networks of deep learning [33].
The selection of training samples for impervious surfaces is mainly based on manual participation in the labor-intensive and time-consuming task, which requires strong prior knowledge and cannot be efficiently applied to automatic impervious surface extraction [34], [35]. Some articles have also trained samples in a semiautomatic manner by computing-related impervious surface indices or integrating ancillary data products while demarcating training samples manually [36], [37], [38]. Additionally, some articles have focused on threshold detection, saliency analysis, and building indices from various remote sensing images for automatic sample extraction, but their practicalities are relatively insufficient and are constrained by the features of objects and pixels, and indices of impervious surfaces. Wang and Li [39] proposed an urban impervious surface automatic threshold detection model from multitemporal Landsat images to reduce manual interference for threshold fine-tuning and reselection, which depends on the capability of the urban impervious surface index to express or depict different land cover types. Tan et al. [40] used two groups of threshold ranges instead of one threshold to obtain more evenly distributed training samples with the parameter values settled empirically through trial and error. Hou et al. [41] designed a saliency-guided edge-preservation framework for semiautomatically updating building databases with limited manual annotation. Zhou et al. [42] developed an automatic sample extraction strategy for both no-change and change samples; however, the strategy is not suitable for impervious surfaces containing high-rise buildings, and it is difficult to determine the parameters for samples with different feature types. Sebari and He [43] proposed an automatic determination of the brightness index threshold based on the histogram frequencies of images to define the discriminative thresholds, but it may confuse some objects presenting the same spectral responses. Liu et al. [44] proposed an automated extraction method of built-up areas to represent the vertical properties of buildings based on multiview images, but it leads to errors from similar features.
Previous articles on remote sensing products of impervious surfaces and land use have mostly used remote sensing images with medium-and-high spatial resolution. Specifically, many products are at medium resolution (10-100 m) with time series, but few are at high resolution (<10 m) with temporal updates [45]. Considering the requirements of refined urban construction and management, the existing numerous remote sensing products, such as the 500-m resolution MODIS data in 2000 [46], 300-m resolution global land cover and land cover change for each year from 1992 to 2015 [47] [55], [56], 10-m resolution European Space Agency (ESA) WorldCover in 2020 [57], and 2-m resolution impervious surface in mainland China in 2017 [58], cannot completely satisfy the needs of spatiotemporal refinement in recent years. The production process of training samples for existing impervious surface products mostly requires manual intervention, which increases the workload and reduces efficiency. In addition, due to the various feature types of ground objects for impervious surfaces, the composition of the urban landscape is complex, especially in areas where different ground objects are mixed. The mixed pixel problem of high-resolution remote sensing images greatly affects the classification results of extracting impervious surfaces, and small ground objects such as houses and roads have a high probability of being wrongly classified [4], [59], [60], [61]. Therefore, an important problem that urgently needs to be solved is how to automatically generating training samples from medium-and-high resolution remote sensing images and products to reduce repetitive manual participation in sample acquisition and improve the efficiency of impervious surface extraction.
To solve the abovementioned problems, this article proposed a novel impervious surface extraction method based on automatically generating training samples from multisource remote sensing products. The goals of this article include the following: 1) propose a method that can automatically generate training samples with several feature types of ground objects; 2) extract impervious surfaces with higher spatial resolution in different years based on a certain limited year and evaluate the classification accuracy; and 3) investigate the feasibility, effectiveness, and advantages of the method proposed for the automatic extraction of impervious surfaces.

A. Study Area
The study area is Wuhan city, which is located in Hubei Province of China with a geographical location of 113°41 E to 115°05 E and 29°58 N to 31°22 N and covers an area of approximately 8569.15 km 2 . It is an economically developed region with a high level of urbanization and rapid urban development. There are various types of interior objects in Wuhan with obvious urban boundaries, which form a concentrated and wide distribution of impervious surfaces. Therefore, automatic extraction of training samples and accurate estimation of impervious surfaces in Wuhan are of vital importance to urban planning, spatial coordination, and environmental protection, but they have rarely been investigated in previous article. In this article, the location of the study area and the spatial distribution of land use types from China's multiperiod LULC remote sensing monitoring dataset (CNLUCC) [62] of Wuhan in 2018 are shown in Fig. 1. The urban land, rural settlement, and other construction land are the typical types of impervious surfaces in CNLUCC.

1) Products and Preprocessing:
The data used in this article mainly include the Sentinel-2 Level-1C product from the ESA [63], global artificial impervious area (GAIA) from Tsinghua University [64], [65], [66], [67], annual urban land maps (AULM) from Liu et al. [68], CNLUCC from the Resource and Environment Science and Data Center of the Chinese Academy of Sciences, and Advanced Land Observing Satellite (ALOS) World 3-D-30 m (AW3D30) version 2.1 product of the global digital surface model (DSM) from the Japan Aerospace Exploration Agency (JAXA) [69]. The seasonal data with strong vegetation growth from June 1 to October 31 are composited to Sentinel-2 image bands by spectral indices to reduce the interference of bare soil for impervious surfaces. The type, spatial resolution, and time coverage of the multisource remote sensing products used in this article are shown in Table I. In addition, the bands and parameters of Sentinel-2 images can be obtained in detail from the Earth Observing System website (https://eos.com/find-satellite/sentinel-2/). The vector data used in this article come from the boundary data of China's counties and cities, including the administrative boundary of Wuhan, and it is mainly used for image cropping.
Sentinel-2 image preprocessing is performed through the Google Earth Engine, including the cloud removal using QA60 band with cloud mask information and the atmospheric correction. For the classification system of products, similar to the GAIA and AULM, which both include two types of impervious surfaces and pervious surfaces, it is necessary to reclassify CNLUCC with six first-level types into two types. The five types of arable land, woodland, grassland, waters, and unused land are reclassified into pervious surfaces, while the type of urban and rural, residential, industrial, and mining land is reclassified into impervious surfaces. The grids representing impervious surfaces and pervious surfaces in the three remote sensing products are assigned values of 1 and 0, respectively, to maintain data consistency.
2) Validation Sample Set: Based on the high-resolution images on the Google Earth platform, the verification sample sets of Sentinel-2 images in Wuhan contain 1671 samples (410 in 2020, which are marked through visual interpretation. To better evaluate the performance of the classification results, the number of impervious surfaces in the validation sample set is relatively large, the numbers of arable land, woodland, and grassland are moderate, the number of easily distinguishable water bodies is small, and unused land with a small total area is the smallest. In addition, due to the large difference in the size and area of a single object for different ground features, all samples are selected to maintain the same size.

III. METHODOLOGY
The proposed impervious surface extraction method based on automatically generating training samples from multisource remote sensing products is presented in Fig. 2. First, a preliminary sample area based on the classification consistency area of remote sensing products and homogeneous area detection of Sentinel-2 images is constructed with impervious surfaces and pervious surfaces. Second, spectral index selection and variance purification calculation are used to improve the data processing efficiency and sample quality after sample overlay analysis. A training sample set containing six types of feature types for ground objects is generated. Third, based on the generated training samples in a certain year with available classification products, the training samples in other years with limited classification products are generated by training sample migration. Finally, the impervious surfaces in a certain year and other years are extracted by random forest model training.

A. Preliminary Sample Area Construction 1) Classification Consistency Area:
To comprehensively consider the consistency and difference of multisource remote sensing products, overlay analysis is used to objectively combine the classification results of GAIA, AULM, and CNLUCC. The classification consistency area is regarded as the true value of the classification result. For the overlay analysis, the original images used by GAIA, AULM, and CNLUCC are all Landsat, and the raster position is the same. By using raster overlay calculation to extract the three products with the same value, we further set the impervious surface grid value to 1 and the pervious surface grid value to 0; that is, the classification consistency area of the existing products can be obtained. Fig. 3 shows the distribution of the classification consistency area of the three products for GAIA, AULM, and CNLUCC of Wuhan in 2018. The red and blue areas in Fig. 3 represent the consistent areas that are classified as impervious surfaces and pervious surfaces for all three products, respectively, which are suitable for scenarios where the types of features are relatively singular and there is a wide distribution of certain features.
By comparing the number of grids with the three areas, more than 80% of the grids have consistent classification results, indicating that the classification results of the three products have relatively small differences. In addition, the yellow area in Fig. 3 indicates the area where the classification results of the three products are inconsistent, which are mostly in areas where populations and buildings are densely populated with other surrounding features. The feature types in the inconsistent classification area are complex and need to be further extracted to improve the classification accuracy of impervious surfaces.
2) Homogeneous Area Detection: To improve the classification accuracy of impervious surfaces, four spectral indices, the normalized difference vegetation index (NDVI) [70], modified normalized difference water index (MNDWI) [71], normalized difference building index (NDBI) [72], and soil extraction index (SoEI) [73], are calculated by the original bands of the Sentinel-2 image. The annual maximum NDVI is used to reflect the surface vegetation density and plant growth. The annual median MNDWI can improve the effectiveness of water extraction, and is beneficial to the case of many buildings in urban areas. The annual minimum NDBI is usually used to extract the urban building land, in particular, to distinguish the spectral characteristics between buildings and bare land. The annual median SoEI is used as the basis for classifying bare soil to distinguish bare soil and impervious surfaces. The elevation of DSM data is also fused to the band of the Sentinel-2 image to provide ground object elevation information of urban buildings and vegetation. Therefore, adding the 13 bands of the original Sentinel-2 image, the composite Sentinel-2 image contains a total of 18 bands.
The ISODATA algorithm is used for the unsupervised classification of the composite Sentinel-2 image, which is divided into n-by-n pixel blocks through cluster analysis. The areas where the pixels in the block of n-by-n pixels are all divided into one type are regarded as homogeneous areas. The spatial resolution of the Sentinel-2 image is 10 m, while the three remote sensing products of GAIA, AULM, and CNLUCC are all 30 m. Therefore, one product grid of GAIA, AULM, or CNLUCC corresponds to 3-by-3 Sentinel-2 image pixels. To avoid mixed pixels with different feature types in the Sentinel-2 image for the sample area in the product grid, it can be judged by calculating whether the 3-by-3 Sentinel-2 image pixels in each product grid belong to the same type. The 3-by-3 Sentinel-2 image pixels contained in the 30 m-by-30 m product grid include two cases (Fig. 4): 1) the nine pixels are of the same type and 2) the nine pixels are not of the same type. For the detection of the homogeneous area for the 10-m resolution Sentinel-2 image with the 30-m resolution remote sensing products, the area that meets the first situation is regarded as a homogenous area that can be ultimately retained; otherwise, it is regarded as a nonhomogeneous area that should be removed.
Homogeneous area detection can be achieved by the Block Statistics method of the spatial analyst in the ArcGIS software [74]. The type of block is set to a 3-by-3 rectangle to calculate the  number of feature types contained in each block in the Sentinel-2 image classification result. The composite Sentinel-2 image of Wuhan in 2018 is calculated by the Block Statistics method, and the pixels with a value of 1 are the homogeneous area to be extracted (Fig. 5). The pixel proportion of the homogenous area is 94.39% of the total pixel number of the composite Sentinel-2 image of Wuhan in 2018, accounting for the majority. Homogeneous areas are generally distributed in the internal area of a single object, especially in areas such as arable land and waters that are widely distributed and occupy a large area, with the pixels forming a large continuous area. In densely populated urban areas and rural areas, the distribution of homogeneous area pixels is significantly more scattered.

1) Preliminary Sample Overlay Analysis:
Through the overlay analysis of the classification consistency area of the remote sensing products and the homogeneous area of the composite  Sentinel-2 image, the intersecting area is the required preliminary sample area. However, the preliminary sample area is still large, and there are some samples with incorrect types that are not consistent with the actual feature types (Fig. 6). Therefore, to further refine the selection of training samples, the preliminary sample area with two types of impervious surfaces and pervious surfaces is analyzed overlapping with the CNLUCC with six feature types.
Through the overlay analysis of the preliminary sample area and the CNLUCC, the type code values of six feature types in CNLUCC are assigned to the corresponding pixels in the same location of the preliminary sample area. The type code values are 1 to 6, corresponding to the first level of CNLUCC (Table II). Among them, the original impervious surface corresponds to CNLUCC type code 5, that is, urban and rural, residential, industrial, and mining land. The original pervious surface corresponds to the other five CNLUCC type codes.
2) Spectral Index Selection: To further improve the accuracy of the training samples, the overlapped pixels in the preliminary sample area are selected and purified by spectral indices. Four spectral indices, NDVI, MNDWI, NDBI, and SoEI, are used to set thresholds to exclude samples with obvious misclassification. According to the field investigation and image comparison, the  Table III. For all image pixels in the preliminary sample area, based on the specific feature type to determine whether the corresponding spectral index meets the filtering threshold conditions, the pixels that do not meet the conditions are removed from the training sample set.

3) Variance Purification Calculation:
After the filtering of spectral index selection, the number of remaining samples is still large and needs to be further purified. The method of variance purification [75], [76], [77] can remove pixels with large deviations from the average value and eliminate incorrect samples. The calculation process of variance purification was as follows: Step 1: Build the average spectral vector (x i ), i = 1, 2, 3, . . . , n of the training samples for each feature type in the study area, x i is the spectral vector, and n is the size of the training samples built from spectral index selection.
Step 2: Calculate the variance Var(x i ) between the spectral vector x i by each pixel for each feature type and the corresponding average spectral vector (x i ), as in Step 3: Calculate the average of the variances aveVar(x i ).
Step 4: Judge each pixel for each feature type to remove impure pixels. If the variance Var(x i ) of the spectral vector of a pixel is greater than A times the average of the variances aveVar(x i ), the pixel can be considered impure. The purified samples p i are calculated as in (2). A is the coefficient that depends on the specific situation to ensure sufficient samples 2, 3, . . . , n. (2) Step 5: Extract the purified training sample set of {p i }, i = 1, 2, 3, . . . , m, where m is the size of purified training samples.
For the average spectral vector, four spectral indices and DSM elevation are merged into the bands of the sentinel-2 image. It is necessary to normalize the samples with all 18 bands of the composite Sentinel-2 image and calculate the average spectral vector and variance. Considering the characteristics of ground objects, different coefficients A of the six feature types are set separately. The urban and rural, residential, industrial, and mining land for impervious surfaces is often easily confused with arable land. Additionally, given the training sample sizes of feature types, the coefficients A of the urban and rural, residential, industrial, and mining land and the arable land are set to 3.5 and 3, respectively, to improve their classification accuracy. The feature types of woodland, grassland, waters, and unused land have obvious characteristics with small gaps, and their coefficients A are all set to 1.5. The training sample set containing six feature types of Wuhan in 2018 is obtained with automated extraction to reduce repetitive manual participation. The quantity and proportion of the six feature types for training samples before and after sample selection and purification are shown in Table IV.

C. Training Sample Migration
Due to the limited coverage time of the remote sensing products of GAIA, AULM, and CNLUCC, only the training samples of Wuhan in 2018 are obtained. However, considering that the actual features may change in different years, the training samples in 2018 cannot be directly used to extract the impervious surfaces in other years. Therefore, to use the constructed training sample set in 2018 for impervious surface extraction over a longer period, the training sample migration method is adopted to extract the training sample sets of Sentinel-2 images of Wuhan in 2019 and 2020. For the composite Sentinel-2 images of Wuhan in 2019 and 2020, the change detection method based on principle component analysis is used to retain the unchanged areas in the sample area. Through overlay analysis, the unchanged areas in 2019 and 2020 are superimposed with the sample area in 2018 to generate the training sample set for the corresponding year.
The specific process of training sample migration is shown in Fig. 2(c). For the other years with limited classification products of impervious pervious or land use (e.g., Wuhan in 2019 and 2020), the composite Sentinel-2 images of this year and a certain year with available classification products (e.g., Wuhan in 2018) are subtracted pixel by pixel and band by band, and the absolute values of the difference images are taken. The main difference information between the difference images is concentrated in  IV  COMPARISON OF THE QUANTITY AND PROPORTION OF THE SIX FEATURE TYPES FOR THE TRAINING SAMPLES OF WUHAN IN 2018 BEFORE AND AFTER SAMPLE  SELECTION

D. Random Forest Model Training
The random forest model has advantages for large and highdimensional data classification; thus, it is selected as the classifier for impervious surface automatic extraction. The input data of random forest model training come from the training samples containing six feature types of Wuhan from 2018 to 2020. Since the input training samples have six feature types, the output classification results also contain the same types. Therefore, the five feature types of arable land, woodland, grassland, waters, and unused land are further reclassified as pervious surfaces, and the one feature type of urban and rural, residential, industrial, and mining land is reclassified as impervious surfaces.
Through the stratified sampling method, 5000 samples of each feature type are selected each year as the training samples of the random forest model. By adding the label of the feature type to each sample, the spectral data of the composite Sentinel-2 images are used as the classification features with a total of 18 features, including the 13 original bands and the added NDVI, MNDWI, NDBI, and SoEI indices and DSM elevation data. The classification features are set as each band of the input image, and the default setting is the square of the total number of features. The feature importance normalization parameters can be returned with total values of 1 through the random forest classifier, which can represent the importance of each feature.
The normalized values of the importance of band features for the composite Sentinel-2 images by random forest model training are shown in Fig. 7. The larger the normalized value is, the greater the importance of the band feature. Among the 18 bands, the B1 (coastal aerosol), B2 (blue), B9 (water vapor), and DSM bands are of higher importance, while the band importance Based on the established validation sample set, the accuracy of the extraction results of the impervious surface distribution in Wuhan from 2018 to 2020 is evaluated (Table V) by the confusion matrix with the overall accuracy (OA), kappa coefficient, user's accuracy (UA), and production's accuracy (PA). With the classification results in this article, the extraction accuracies of impervious surfaces for Wuhan in 2018, 2019, and 2020 were 94.02%, 94.45%, and 93.87%, respectively. This shows that the training samples obtained by the method proposed in this article have high accuracy, and the classification of impervious surfaces performs well. Thus, the training sample sets of impervious surfaces in existing years generated according to the method proposed can be used for Sentinel-2 images or other high-resolution images in other years by training sample migration for long-term series of impervious surface extraction.

2) Comparison With Existing Remote Sensing Products:
Compared with the classification results of the two 10-m resolution remote sensing products of the ESA WorldCover 2020 and ESRI 2020 LULC, and the impervious surfaces of Wuhan in 2020 extracted in this article, the feasibility and effectiveness of the method proposed are discussed. Six local areas where the types of more complex ground objects at the junction of urban and rural areas are selected under the original remote sensing image of Sentinel-2 Level-1C production are shown in Fig. 9. The results show that the degree of refinement of the ESRI 2020 LULC is relatively low, and it has an overestimation of the impervious surfaces in case of the pervious surfaces exist in the impervious surface area. For the ESA WorldCover 2020, the road boundaries of houses and roads with a small area can be accurately extracted with fewer missing points, but it has poor extraction of bright impervious surfaces, resulting in an underestimation of impervious surfaces. For the result in this article, it is more conducive to accurately distinguishing the boundary between impervious and pervious surfaces, the relatively obvious boundaries of houses and roads, and it also has a more refined extraction of bright impervious surfaces. The accuracy assessments of the impervious surface extraction results obtained by this article, ESA WorldCover 2020, and the ESRI 2020 LULC are shown in Fig. 10. Therefore, the impervious surface training samples extracted by the proposed method are practical. B. Discussion 1) Accuracy Assessment of the Generated Training Samples: In terms of sample accuracy assessment, Olofsson et al. [78] pointed out that probability sampling design can be implemented considering actual conditions. Given that the number of sample sets is still large after sample selection and purification, this article randomly selects the same number of each feature type to improve efficiency. A total of 5000 samples of each feature type are randomly selected by the stratified random sampling method. Additionally, 5% of the stratified random samples, that  is, 250 samples of each feature type, are also randomly selected for accuracy assessment. Through visual interpretation of Sentinel-2 images, the confusion matrix [79] with OA, kappa coefficient, UA, and PA is used to evaluate the accuracy of the classification results for the training samples. The accuracy assessment results of the six feature types for the stratified random samples of Wuhan in 2018 are displayed in the confusion matrix with an OA of 91.07% (Table VI).
To evaluate the sample accuracy of the impervious surface and pervious surface, the six feature types are reclassified into two types. Among them, five feature types of arable land, woodland, grassland, waters, and unused land are regarded as pervious surfaces, and the feature type of urban and rural, residential, industrial, and mining land is regarded as impervious surfaces. The accuracy assessment results of the two types of stratified random samples of Wuhan in 2018 are displayed in the confusion matrix with an OA of 99.4% (Table VII). Therefore, it is conducive to the performance of the training samples for impervious surfaces in practical applications.
2) Sample Separability Assessment: To verify the separability of training samples for the impervious surfaces and the background pervious surfaces, both the Jeffries-Matusita (JM) distance [80] and the transformed divergence (TD) distance [81] are calculated in this article. The JM distance and the TD distance for the training samples based on remote sensing data of Wuhan in 2018 are 1.84 and 1.99, respectively. The JM distance of impervious surfaces is more than 1.8, which indicates the high separability between the impervious surfaces  VII  ACCURACY ASSESSMENT WITH THE CONFUSION MATRIX OF TWO TYPES, I.E., IMPERVIOUS SURFACES AND PERVIOUS SURFACES, FOR THE STRATIFIED RANDOM  SAMPLES FOR WUHAN IN 2018 and the background pervious surfaces, and the powerful capability of the automatically generating training sample's method for impervious surface extraction.
3) Limitations and Future Article: In this article, Sentinel-2 images are directly used as the original data, and multiple spectral indices, including NDVI, MNDWI, NDBI, and SoEI, as well as DSM elevation data, are added as auxiliary data. Only the spectral information of the composite Sentinel-2 images is used in the classification process. Next article can consider adding other features, for example, texture, to improve the efficiency and accuracy of impervious surface classification. Meanwhile, to further improve the accuracy of impervious surfaces, the next article will also consider the distinction between impervious surface types, including the bright impervious surfaces and dark impervious surfaces. In addition, the experiment in this article only conducts a case study on Wuhan city in 3 years, which verifies the feasibility of the proposed method. In the future, it will be explored in adopting other high-resolution images with longer time coverage for the long-term series analysis of impervious surfaces and improving the applicability of spectral index thresholds in different regions.

V. CONCLUSION
It is important to reduce manual participation in training sample acquisition and improve the efficiency of impervious surface extraction. The purpose of this article is to propose a method of impervious surface extraction based on automatically generating training samples from Sentinel-2 images, multisource remote sensing products of GAIA, AULM, and CNLUCC, as well as DSM data. The contributions of this article are primarily twofold: 1) a method is proposed for automatically generating training samples of impervious surfaces with 10-m resolution and improved sample quality, and 2) a powerful strategy is designed through training sample migration, which can be applied to obtain the training samples in the year with limited classification products based on the certain year with available classification products. The accuracy assessment of the training samples showed that the extracted training sample set of Wuhan in 2018 had a higher accuracy with an OA of 99.40% and kappa coefficient of 0.98 with the automatic extraction process, thus reducing repetitive manual participation. In addition, due to the time constraints of remote sensing products for impervious surfaces and land use, based on the generated training sample set, sample migration was adopted to perform change detection based on principle component analysis to obtain the respective training sample sets in 2019 and 2020. To further validate the availability of the training samples by the method proposed, the obtained training sample sets were used for random forest model training to extract the distribution of impervious surfaces of Wuhan in 2018, 2019, and 2020 with OAs of 94.02%, 94.45%, and 93.87%, respectively. The spatial resolution of the final classification result of the impervious surfaces in this article was improved from 30 m of the original remote sensing products to 10 m. Yuqin Wu received the B.Sc. degree in geographic