Multi-Level Classification Based on Trajectory Features of Time Series for Monitoring Impervious Surface Expansions

As urbanization has profound effects on global environmental changes, quick and accurate monitoring of the dynamic changes in impervious surfaces is of great significance for environmental protection. The increased spatiotemporal resolution of imagery makes it possible to construct time series to obtain long-time-period and high-accuracy information about impervious surface expansion. In this study, a three-step monitoring method based on time series trajectory segmentation was developed to extract impervious surface expansion using Landsat time series and was applied to the Xinbei District, Changzhou, China, from 2005 to 2017. Firstly, the original time series was segmented and fitted to remove the noise caused by clouds, shadows, and interannual differences, leaving only the trend information. Secondly, the time series trajectory features of impervious surface expansion were described using three phases and four types with nine parameters by analyzing the trajectory characteristics. Thirdly, a multi-level classification method was used to determine the scope of impervious surface expansion, and the expansion time was superimposed to obtain a spatiotemporal distribution map. The proposed method yielded an overall accuracy of 90.58% and a Kappa coefficient of 0.90, demonstrating that Landsat time series remote sensing images could be used effectively in this approach to monitor the spatiotemporal expansion of impervious surfaces.


Introduction
Over the past few decades, human activities have changed the global environment at an unprecedented rate and scale [1].For example, urbanization has rapidly increased the impervious surface area and has degraded large swaths of high-quality farmland and ecological settings [2][3][4].As this impervious surface expansion profoundly affects food security and the ecosystems around the world [5,6], monitoring this trend is important for global sustainable development.The rapid development of Earth observation technology has made remote sensing data increasingly abundant and applicable to the monitoring of impervious surface expansion [7].
Many methods have been developed to monitor impervious surface expansion in recent years, such as impervious surface index analysis, spectral mixture analysis, and random forest analysis [8][9][10].
The corresponding studies have been focused on finding spectral differences between impervious surfaces and other land covers using single or multiple images.However, due to the presence of mixed pixels in coarse or medium-resolution images, relying solely on multi-temporal data often does not yield satisfactory results [11,12].Time series images can better solve the spectral confusion problem between similar land cover types [13][14][15][16] and provide clearer distinctions between vegetation changes caused by phenological changes, interannual variations, and vegetation changes caused by impervious surface changes [17,18].Time series change detection through two main approaches has become increasingly popular [19,20]: time series similarity measurement and trajectory classification [21,22].In the former method, the similarity between the time series of target pixels and a reference time series consisting of sample pixels of various land cover types is calculated, and the land cover type of the sample pixels most similar is set to the target pixels [23].This method requires higher quality time series for the reference samples and makes it necessary to select sample points, including all times and change types [24].Therefore, sample selection is difficult or even impossible for areas in which land use changes frequently over time.Meanwhile, in trajectory classification, the features of different types of land cover changes are extracted from time series trajectories, and the time series trajectories are classified based on these features using classification methods [25][26][27][28].Trajectory classification reduces the difficulty caused by different change times in similarity measurements.Two types of trajectory classification algorithms, hypothesized trajectory and multi-date classification algorithms, have many types of applications and yield good results [22].
Although trajectory classification provides a promising means of analyzing impervious surface expansion, several challenges remain.First, noises (such as clouds, shadows, and interannual differences) make the original time series extremely complicated and obscure the characteristics of the expansion trajectory.It is difficult to extract valuable trajectory features.Second, the time series trajectory of impervious surface expansion is complex and difficult to depict using parameters, which makes it necessary to screen the feature parameters and establish a remote sensing information extraction model.
Considering the shortcomings mentioned above, we proposed a multi-level classification approach combining support vector machine (SVM) and decision tree classification based on the trajectory features of the time series.We generated the time series trajectory of each pixel and extracted the features of impervious surface expansion trajectories, then employed multi-level classification to obtain impervious surface expansion information.To verify the performance, the proposed approach was utilized to detect impervious surface expansion in Xinbei, Changzhou from 2005 to 2017.

Study Area and Data
Xinbei District in northern Changzhou, Jiangsu, China (31 • 48 -32 • 03 N; 119 • 46 -120 • 01 E; Figure 1) was selected as the study area.Situated in the Yangtze River Delta, this area was among the first national high-tech industrial development zones in China (established in 1992) and has developed rapidly due to its flat terrain, convenient transportation access, and superior geographical conditions.Xinbei District has jurisdiction over three subdistricts (Hehai, Sanjing, and Longhutang) and six towns (Chunjiang, Menghe, Xixiashu, Luoxi, Xuejia, and Xinqiao), and its annual GDP growth rate has exceeded 20% for many years.In 2006, the Changzhou Municipal Government moved to Xinbei District, making this district one of the key development areas in Changzhou.In the past 10 years, there has been great demand for development land in Xinbei District, leading to numerous expansion projects.The remote sensing data used in this study were downloaded from the United States Geological Survey (USGS) website and the Geospatial Data Cloud.The image strip number was 119, and the row number was 38.Three types of satellite remote sensing data from 2005 to 2017 were collected: Landsat5 Thematic Mapper (TM), Landsat7 Enhance Thematic Mapper Plus (ETM+), and Landsat8 Operational Land Imager (OLI).To avoid detection errors caused by inconsistent crop growth periods, we selected images acquired during periods of thriving crop growth every year (Table 1).The land use status data were obtained from the Changzhou Land and Resources Bureau, including land use change survey data in Changzhou City from 2005 to 2016.These data included 10 types: cultivated land, garden areas, woodland, grassland, water, transport land, urban areas, industrial land, other agricultural land, and scenic land.The impervious surface includes transport land, urban land, and industrial land.The non-impervious surface includes cultivated land, garden, woodland, grassland, water, other agricultural land, and scenic land.The remote sensing data used in this study were downloaded from the United States Geological Survey (USGS) website and the Geospatial Data Cloud.The image strip number was 119, and the row number was 38.Three types of satellite remote sensing data from 2005 to 2017 were collected: Landsat5 Thematic Mapper (TM), Landsat7 Enhance Thematic Mapper Plus (ETM+), and Landsat8 Operational Land Imager (OLI).To avoid detection errors caused by inconsistent crop growth periods, we selected images acquired during periods of thriving crop growth every year (Table 1).The land use status data were obtained from the Changzhou Land and Resources Bureau, including land use change survey data in Changzhou City from 2005 to 2016.These data included 10 types: cultivated land, garden areas, woodland, grassland, water, transport land, urban areas, industrial land, other agricultural land, and scenic land.The impervious surface includes transport land, urban land, and industrial land.The non-impervious surface includes cultivated land, garden, woodland, grassland, water, other agricultural land, and scenic land.

Overview of Methods
Based on the Landsat time series fitting trajectory, a multi-level classification method involving a combination of SVM and decision tree classification was proposed, and an impervious surface expansion monitoring framework was established in which time series segmentation was fused with SVM and decision tree classification.The flowchart of this study is shown in Figure 2. Firstly, the time series segmentation algorithm was used to obtain the trajectory of the spectral changes of each pixel.Next, the fitting trajectories of the sample points were analyzed to extract the features of the time series trajectories of impervious surface expansion.Subsequently, the multi-level classification method was used to extract the scope of impervious surface expansion.Finally, the scope and time of impervious surface expansion were superimposed to obtain a spatiotemporal distribution map of impervious surface expansion.Based on the Landsat time series fitting trajectory, a multi-level classification method involving a combination of SVM and decision tree classification was proposed, and an impervious surface expansion monitoring framework was established in which time series segmentation was fused with SVM and decision tree classification.The flowchart of this study is shown in Figure 2. Firstly, the time series segmentation algorithm was used to obtain the trajectory of the spectral changes of each pixel.Next, the fitting trajectories of the sample points were analyzed to extract the features of the time series trajectories of impervious surface expansion.Subsequently, the multi-level classification method was used to extract the scope of impervious surface expansion.Finally, the scope and time of impervious surface expansion were superimposed to obtain a spatiotemporal distribution map of impervious surface expansion.

Preprocessing of Landsat Images
Geometric registration was completed for Landsat Level-1 data products obtained from the Geospatial Data Cloud and USGS website, with errors within 0.5 pixels, which was sufficiently accurate for this study [29].The Quick Atmospheric Correction method in the Environment for Visualizing Images (ENVI, a software for remote sensing image processing) was used to perform

Preprocessing of Landsat Images
Geometric registration was completed for Landsat Level-1 data products obtained from the Geospatial Data Cloud and USGS website, with errors within 0.5 pixels, which was sufficiently accurate for this study [29].The Quick Atmospheric Correction method in the Environment for Visualizing Images (ENVI, a software for remote sensing image processing) was used to perform atmospheric corrections for the Landsat TM/ETM+/OLI data.Due to the Landsat7 server failure in June 2003, the images acquired later were stripped, and the strips were recovered by replacing them with available pixels at adjacent time points [30,31].For the connectivity and integrity of feature boundaries of remote sensing images, a k-nearest neighbor (KNN) filter was used to filter each image [32,33].The images processed by atmospheric correction, strip recovery, and KNN filtering were clipped to obtain images of the study area.In addition, the normalized difference vegetation index (NDVI) was calculated for each of the images processed as described above, and NDVI time series were then constructed to prepare for impervious surface expansion monitoring.

Trajectory Generation
Although the remote sensing data sources were screened to obtain images of the same season each year, the time series constructed had extremely high and low peaks due to the presence of a few clouds and shadows in the original images; slight interannual differences were also found in the time series data due to different imaging conditions, such as light and moisture in the atmosphere for each image.These factors all reduced the impervious surface expansion monitoring accuracy and were collectively referred to as time series noise.Although the filter-based noise reduction method has been greatly improved in recent years, the problems of missing data and excessive smoothing in time series have not been resolved [34,35].To eliminate such noise and obtain the spectral changes of each pixel, the LandTrendr algorithm was employed, which has provided accurate results when employed to detect forest disturbances, farming abandonment in agricultural land, and disturbance-recovery of open-pit mines [36][37][38][39].
As shown in Figure 3, the LandTrendr algorithm included a segmentation strategy before fitting to eliminate noise while preserving as much of the change information of the original time series as possible.Firstly, the entire curve was iterated multiple times by minimizing the residuals, and the positions of the segmentation points in the time series trajectory were identified.To prevent iteration from causing over-fitting, the angle threshold was further discriminated after segmentation.Thereby, the excess segmentation points were removed, and the remaining vertices were used as the segmentation vertices.Then, based on the principle of the mean-square error, the segmentation vertices were linearly fitted to obtain a group of connected segments as the time series trajectory of each pixel.

Feature Extraction of Impervious Surface Expansion Trajectories
The process of impervious surface expansion changes with time.The expansion may occur at any time point, but the features of the time series trajectories are consistent before and after the expansion time point.In this study, the time dimension was reduced.The segments far from the change time point were ignored, and only the short time trajectories before, during, and after the change were considered.The features of these short time trajectories were analyzed, and the template for the impervious surface expansion trajectories was developed based on the results.The classification algorithm was used to find the pixels matching this template and to obtain the scope of impervious surface expansion.
As shown in Figure 4, the fitting trajectories of the impervious surface expansion time series were represented by three phases, four types, and nine parameters.Most of the pixels corresponding to impervious surface expansion first underwent a relatively stable period in which the NDVI fluctuated around L1 (I).After this period, the NDVI dropped substantially (II), then continued to decline or rose slightly.However, the magnitude of this decline or recovery was small.Finally, the NDVI became steady at L2 (III).A few pixels underwent only one or two phases, but their NDVIs fell from a relatively undetermined level (L1) to a stable level (L2) as well, as shown in Figure 4b-d.

Feature Extraction of Impervious Surface Expansion Trajectories
The process of impervious surface expansion changes with time.The expansion may occur at any time point, but the features of the time series trajectories are consistent before and after the expansion time point.In this study, the time dimension was reduced.The segments far from the change time point were ignored, and only the short time trajectories before, during, and after the change were considered.The features of these short time trajectories were analyzed, and the template for the impervious surface expansion trajectories was developed based on the results.The classification algorithm was used to find the pixels matching this template and to obtain the scope of impervious surface expansion.
As shown in Figure 4, the fitting trajectories of the impervious surface expansion time series were represented by three phases, four types, and nine parameters.Most of the pixels corresponding to impervious surface expansion first underwent a relatively stable period in which the NDVI fluctuated around L1 (I).After this period, the NDVI dropped substantially (II), then continued to decline or rose slightly.However, the magnitude of this decline or recovery was small.Finally, the NDVI became steady at L2 (III).A few pixels underwent only one or two phases, but their NDVIs fell from a relatively undetermined level (L1) to a stable level (L2) as well, as shown in Figure 4b-d.To consider more trajectory features, we took the impervious surface expansion time series trajectory undergoing all three phases as an example.The features of the other three types of trajectories could be regarded as subsets.The features of the time series trajectory of impervious surface expansion were expressed using the parameters listed in Table 2.

Multi-Level Classification for Time Series Trajectories
There is a shape difference between the trend trajectory of the impervious surface expansion pixels and those of the other land use types.With this difference, the classification method can be used to identify the impervious surface expansion pixels.Therefore, it is necessary to select an effective classification method to classify time series trajectories.To consider more trajectory features, we took the impervious surface expansion time series trajectory undergoing all three phases as an example.The features of the other three types of trajectories could be regarded as subsets.The features of the time series trajectory of impervious surface expansion were expressed using the parameters listed in Table 2.

Multi-Level Classification for Time Series Trajectories
There is a shape difference between the trend trajectory of the impervious surface expansion pixels and those of the other land use types.With this difference, the classification method can be used to identify the impervious surface expansion pixels.Therefore, it is necessary to select an effective classification method to classify time series trajectories.
Firstly, we chose to use an SVM for trajectory classification for the following reasons.In recent years, SVM supervised classification has been widely used for remote sensing image classification [40].The basic idea of SVM classification is to transform the actual problem into a high-dimensional feature space and find the best-fit plane in this space, which serves as the classification plane to distinguish different types.It has yielded good classification accuracy in many impervious surface monitoring studies [41].Therefore, we applied SVM classification to identify the impervious surface expansion trajectories.
However, there are various types of land use changes.For example, the time series trajectory of cultivated land changing to pond water will exhibit three phases: steady, declining, and re-steady periods, which is similar to the time series trajectory of impervious surface expansion shown in Figure 4a.However, they are different types of change.Consequently, more accurate impervious surface expansion information is required to screen the pixel trajectories further.SVM is more applicable to binary classification, while the screening process in this case involves four different trajectory types, which is considered to be a multi-class classification issue [41].A decision tree can be used to distinguish many different land types by analyzing the differences between the spectral curves of typical features [42].Schneider showed that decision trees are powerful tools for monitoring land cover changes in suburban areas [11].Moreover, the rules of decision trees are relatively flexible and can be improved based on the classification results.So, the combination of an SVM and a decision tree can improve the classification accuracy [43].
Therefore, SVM and decision tree classification were combined to construct a multi-level classification method to classify the time series trajectories and identify the scope of impervious surface expansion in this study.The multi-level classification process for time series trajectories is shown in Figure 5.The maximum recovery value, maximum disturbance value, and ratio of maximum disturbance to maximum recovery were selected from the trajectory parameters as feature space 1. Impervious surface expansion samples were obtained based on the land use change survey data from 2005 and 2016.SVM classification was employed to extract the pixels corresponding to impervious surface expansion trajectories.The NDVI upon completion of expansion, NDVI upon completion of recovery, and average NDVI after expansion were selected from the remaining trajectory parameters to form feature space 2. With reference to Google Earth high-resolution images, non-impervious surface expansion samples were selected from the SVM extraction results.The box plot method was used to determine the decision tree thresholds, and the decision tree rules were then employed to filter the scope of impervious surface expansion further.Finally, the expansion scope and time were superimposed to obtain the spatiotemporal distribution of impervious surface expansion.Firstly, we chose to use an SVM for trajectory classification for the following reasons.In recent years, SVM supervised classification has been widely used for remote sensing image classification [40].The basic idea of SVM classification is to transform the actual problem into a high-dimensional feature space and find the best-fit plane in this space, which serves as the classification plane to distinguish different types.It has yielded good classification accuracy in many impervious surface monitoring studies [41].Therefore, we applied SVM classification to identify the impervious surface expansion trajectories.
However, there are various types of land use changes.For example, the time series trajectory of cultivated land changing to pond water will exhibit three phases: steady, declining, and re-steady periods, which is similar to the time series trajectory of impervious surface expansion shown in Figure 4a.However, they are different types of change.Consequently, more accurate impervious surface expansion information is required to screen the pixel trajectories further.SVM is more applicable to binary classification, while the screening process in this case involves four different trajectory types, which is considered to be a multi-class classification issue [41].A decision tree can be used to distinguish many different land types by analyzing the differences between the spectral curves of typical features [42].Schneider showed that decision trees are powerful tools for monitoring land cover changes in suburban areas [11].Moreover, the rules of decision trees are relatively flexible and can be improved based on the classification results.So, the combination of an SVM and a decision tree can improve the classification accuracy [43].
Therefore, SVM and decision tree classification were combined to construct a multi-level classification method to classify the time series trajectories and identify the scope of impervious surface expansion in this study.The multi-level classification process for time series trajectories is shown in Figure 5.The maximum recovery value, maximum disturbance value, and ratio of maximum disturbance to maximum recovery were selected from the trajectory parameters as feature space 1. Impervious surface expansion samples were obtained based on the land use change survey data from 2005 and 2016.SVM classification was employed to extract the pixels corresponding to impervious surface expansion trajectories.The NDVI upon completion of expansion, NDVI upon completion of recovery, and average NDVI after expansion were selected from the remaining trajectory parameters to form feature space 2. With reference to Google Earth high-resolution images, non-impervious surface expansion samples were selected from the SVM extraction results.The box plot method was used to determine the decision tree thresholds, and the decision tree rules were then employed to filter the scope of impervious surface expansion further.Finally, the expansion scope and time were superimposed to obtain the spatiotemporal distribution of impervious surface expansion.

SVM for Extracting Pixels with Features of Impervious Surface Expansion Trajectories
For feature space 1, typical areas witnessing impervious surface expansion and no impervious surface expansion were selected as the areas of interest.SVM classification was used to extract the pixels that had features of impervious surface expansion trajectories, which was considered the preliminary extraction of impervious surface expansion.To identify the pixels that had features of impervious surface expansion trajectories, we selected three indices (maximum disturbance value, maximum recovery value, and ratio of maximum disturbance to maximum recovery) as SVM classification features (feature space 1).These three indices were selected because they can basically be used to determine the trend features of trajectories.
The survey data of land use changes from 2005 and 2016 were superimposed to obtain the scope of planned impervious surface expansion.Within this scope, the pixels having the fitting trajectories of the impervious surface expansion time series were uniformly selected as the first category of sample points, and those without any significant changes in the time series were selected as the second category of sample points.The sample point distribution is shown in Figure 6.The sample points in the first category were of four types that corresponded to the four different types of impervious surface expansion trajectory features in Figure 4; those in the second category were pixels that did not have impervious surface expansion, whose corresponding fitting trajectories were horizontal or nearly horizontal lines.Finally, all of the pixels in the study area were designated as pixels with or without the features of impervious surface expansion trajectories.The pixels without the features of impervious surface expansion trajectories were discarded, while those with these features were classified further.
Remote Sens. 2019, 11, x FOR PEER REVIEW 9 of 23 For feature space 1, typical areas witnessing impervious surface expansion and no impervious surface expansion were selected as the areas of interest.SVM classification was used to extract the pixels that had features of impervious surface expansion trajectories, which was considered the preliminary extraction of impervious surface expansion.
To identify the pixels that had features of impervious surface expansion trajectories, we selected three indices (maximum disturbance value, maximum recovery value, and ratio of maximum disturbance to maximum recovery) as SVM classification features (feature space 1).These three indices were selected because they can basically be used to determine the trend features of trajectories.
The survey data of land use changes from 2005 and 2016 were superimposed to obtain the scope of planned impervious surface expansion.Within this scope, the pixels having the fitting trajectories of the impervious surface expansion time series were uniformly selected as the first category of sample points, and those without any significant changes in the time series were selected as the second category of sample points.The sample point distribution is shown in Figure 6.The sample points in the first category were of four types that corresponded to the four different types of impervious surface expansion trajectory features in Figure 4; those in the second category were pixels that did not have impervious surface expansion, whose corresponding fitting trajectories were horizontal or nearly horizontal lines.Finally, all of the pixels in the study area were designated as pixels with or without the features of impervious surface expansion trajectories.The pixels without the features of impervious surface expansion trajectories were discarded, while those with these features were classified further.

Decision Tree for Screening Impervious Surface Expansion Pixels
The SVM extracted curves with the shapes of the impervious surface expansion trajectories.These curves represented the land use changes, but not all of them indicated impervious surface expansion.Therefore, it was necessary to classify the SVM results further.After analyzing the differences between the trajectories, we chose a decision tree classification method to refine impervious surface expansion pixels.
It can be seen from the SVM extraction results that there were two main types of pixels having the shapes of the impervious surface expansion trajectories in the study area: vegetation changing to impervious surface and vegetation changing to pond water.As shown in Figure 4, Phases I and III represent the spectral features of land before and after the change, respectively.The spectral features of the land after the change were analyzed, and the impervious surface expansion pixels were selected from all of the pixels corresponding to impervious surface expansion trajectories.The samples of vegetation changing to impervious surface and of vegetation changing to pond water were selected from the SVM classification results.Then, using a box plot to analyze the statistical characteristics of sample points, the decision tree parameters and the threshold were determined according to the statistical characteristics in the box plot.To obtain the divisibility of the time series trajectory parameters of the samples, we calculated the reparability of the NDVI upon completion of expansion, NDVI upon completion of recovery, and average NDVI after expansion; then, the parameters with better divisibility were selected to construct a decision tree.
In this study, the impervious surface expansion trajectories had four different shapes.Since the parameters for each shape may not have been stored in a single file, it was necessary to find the corresponding parameter file according to the features of the different shapes.The parameters for each segment after applying LandTrendr, such as vertices, mean values, and disturbance values, were stored in the corresponding band.The number of segments and NDVI disturbance value corresponding to each segment were used to find the parameters with better divisibility corresponding to each shape, which were then subjected to decision tree identification.The decision tree constructed for the NDVI upon completion of expansion is shown in Figure 7.After the scope of impervious surface expansion was determined using the method described above, it was also necessary to superimpose the results with the impervious surface expansion time to obtain the spatiotemporal distribution of impervious surface expansion.The calculation of the impervious surface expansion time (represented by parameter i) is described in Section 2.4.Firstly, the maximum disturbance segment of the time series trajectory was detected, and the year following After the scope of impervious surface expansion was determined using the method described above, it was also necessary to superimpose the results with the impervious surface expansion time to obtain the spatiotemporal distribution of impervious surface expansion.The calculation of the impervious surface expansion time (represented by parameter i) is described in Section 2.4.Firstly, the maximum disturbance segment of the time series trajectory was detected, and the year following the disturbance was used as the impervious surface expansion time.This method was employed because the land use type at the beginning of NDVI disturbance was still the type before the change, and the land use type changed formally in the year following the NDVI disturbance.

Impervious Surface Expansion Identification Results
Based on Google Earth high-resolution remote sensing images, 255 samples of vegetation changing to impervious surface and 248 samples of vegetation changing to pond water were selected based on the SVM classification results.Then, the divisibility of the NDVI upon completion of expansion, NDVI upon completion of recovery, and average NDVI after expansion were collected for these samples.
It can be seen from Figure 8 that the NDVI upon completion of recovery and the average NDVI after expansion have poor divisibility and large overlap areas.The NDVI upon completion of expansion has the best divisibility for the two land types, and pond water could be extracted from the SVM results by setting a threshold T. After conducting multiple tests, T was set to 150.If the NDVI upon completion of expansion was greater than 150, the land type was impervious surface that needed to be preserved; if it was less than 150, the land type was non-impervious surface that needed to be filtered out.upon completion of expansion was greater than 150, the land type was impervious surface that needed to be preserved; if it was less than 150, the land type was non-impervious surface that needed to be filtered out.SVM classification was performed to classify the three parameters (NDVI disturbance value, NDVI recovery value, and the ratio of maximum disturbance to maximum recovery) to obtain the pixels with the features of impervious surface expansion trajectories.As shown in Figure 9a, many areas that did not change to impervious surface were identified as exhibiting impervious surface expansion.This situation was quite severe in the northern and western parts of the study area, SVM classification was performed to classify the three parameters (NDVI disturbance value, NDVI recovery value, and the ratio of maximum disturbance to maximum recovery) to obtain the pixels with the features of impervious surface expansion trajectories.As shown in Figure 9a, many areas that did not change to impervious surface were identified as exhibiting impervious surface expansion.This situation was quite severe in the northern and western parts of the study area, specifically, Regions I and II in Figure 9.It can be seen from the high-resolution images that the case of vegetation cover area, such as cultivated land and gardens, changing to pond water was the most confused with actual impervious surface expansion.Because water absorbs weakly in the green band and strongly in the near-infrared band, the NDVI of pond water is low.If vegetation changes to pond water, the NDVI will also undergo three phases: steady, declining, and re-steady periods.The SVM failed to distinguish such pixels.Therefore, a decision tree was used to classify and filter the pixels of non-impervious surface expansion further, with the classification results shown in Figure 9b.It can be seen that the incorrectly classified pixels in Regions I and II were almost completely removed, but there were no significant effects on Regions III and IV, which witnessed actual impervious surface expansion.The decision tree also had obvious removal effects on other incorrectly classified pixels scattered throughout the study area.No distinction was made between bare soil and impervious surface in this study, mainly because Changzhou is now in a period of rapid economic development, and its impervious surface area is expanding on a large scale.Occasionally, small parcels of bare soil are maintained for some time due to the approval process.However, once the process is completed, these areas quickly change into impervious surface.

Evaluation of Impervious Surface Expansion Monitoring Accuracy
Verification points were selected to calculate the confusion matrix to evaluate the accuracy of No distinction was made between bare soil and impervious surface in this study, mainly because Changzhou is now in a period of rapid economic development, and its impervious surface area is expanding on a large scale.Occasionally, small parcels of bare soil are maintained for some time due to the approval process.However, once the process is completed, these areas quickly change into impervious surface.

Evaluation of Impervious Surface Expansion Monitoring Accuracy
Verification points were selected to calculate the confusion matrix to evaluate the accuracy of the impervious surface expansion monitoring results.It can be seen from the final results that the images of the study area were divided into 13 layers, including 12 layers of impervious surface expansion and one layer of no-impervious surface expansion.Using hierarchical random sampling, we randomly selected 520 samples (pixels of 40 samples in each layer) for the accuracy test.The samples were verified by visual interpretation of the 20 issues of Landsat images used in the study, assisted by high-resolution Google Earth images.
It can be seen from Table 3 that the impervious surface expansion monitoring results calculated by using the confusion matrix had an overall accuracy of 90.58% and a Kappa coefficient of 0.90.User accuracy is between 85% and 95%, and producer accuracy is between 85.4% and 94.7%, which suggests that our maps provide good estimates of the spatiotemporal variation of the impervious surface.Thus, the multi-level impervious surface expansion identification method based on time series segmentation could be used to monitor the impervious surface expansion accurately.However, among the 40 sample points with no impervious surface expansion that were found in the monitoring results, five samples witnessed impervious surface expansion in different years, and among the sample points with impervious surface expansion, only three witnessed no impervious surface expansion.Therefore, the impervious surface scope was underestimated to some extent.Overall accuracy: 90.58% Kappa coefficient: 0.90 Total number of reference pixels: 520

Reference
Note: NC denotes no impervious surface expansion.The fill color of yellow denotes samples whose change types were correctly detected.The pink denotes samples whose change types were not detected.
In addition, the yellow area in Table 3 shows the monitoring accuracy of the impervious surface expansion time.The time monitoring error was generally 1 year, but it was 2 years in a few cases.The user and producer accuracies for monitoring the change time are greater than 85% in each year and are greater than 90% in most cases, and the years with less than 90% accuracy are the first and last two to three years in the time series.

Spatiotemporal Features of Impervious Surface Expansion
The spatiotemporal features of impervious surface expansion in Xinbei District are shown in Figure 10.The impervious surface expansion clearly varies by region.Hehai, Sanjing, and Longhutang subdistricts in the south are the economic centers of Xinbei District, and the impervious surface there was near saturation before 2005 and barely expanded after that.Xuejia, Luoxi, and Xinqiao Towns witnessed the widest impervious surface expansion to the periphery, followed by Chunjiang Town in the northeast.The impervious surface expansion scope was the smallest in Menghe and Xixiashu Towns.Besides the regional differences, there was also some regularity in the impervious surface expansion.The urban impervious surface basically presented a circular expansion pattern ranging from the subdistricts and towns as the centers to the periphery.The impervious surface for transportation served as the link between townships and subdistricts, showing a networked expansion pattern.
The impervious surface expansion in Xinbei District also exhibited obvious temporal trends.In 2005-2017, the impervious surface expansion area showed a shrinking trend with alternating peaks and valleys.After 2005, the impervious surface area expanded by 3427.79 hectares, with an average annual expansion of about 270 hectares.Over time, the impervious surface expansion area declined overall, and it fell to a minimum of only 45.54 hectares in 2017.Impervious surface expansion peaks occurred in 2006, 2009, and 2014, with expansion areas of 866.07, 583.65, and 380.34 hectares, respectively, which were significantly higher than those in the adjacent years.

Advantages of Multi-Level Impervious Surface Expansion Monitoring
Multi-level classification was combined with trajectory features of time series to monitor impervious surface expansion.Compared with the existing research, three improvements were made in the proposed approach for detecting the impervious surface expansion time and scope.First, the method of post-classification change detection has higher classification accuracy requirements, and the final change monitoring error is the accumulation of the individual image classification errors [44][45][46].In addition, the time-series-based method avoids manipulation of a single image, thereby reducing the number of intermediate links with possible errors.Second, the time series similarity measurement can be used to identify impervious surface expansion finely, but it relies heavily on the change samples, because it is necessary to provide change samples in different time periods.The proposed approach in this study only requires samples with a certain change feature to determine the type of change at any time point.Third, the method of trajectory classification can indeed be employed to detect impervious surface expansion at different time points based on a small number of samples, but the accuracy is low due to noise because the original time series are used.In this study, noise caused by, e.g., clouds, shadows, and interannual difference, was eliminated, and only the trend trajectory of time series was preserved by segmentation fitting, thereby increasing the impervious surface expansion monitoring accuracy.

Advantages of Multi-Level Impervious Surface Expansion Monitoring
Multi-level classification was combined with trajectory features of time series to monitor impervious surface expansion.Compared with the existing research, three improvements were made in the proposed approach for detecting the impervious surface expansion time and scope.First, the method of post-classification change detection has higher classification accuracy requirements, and the final change monitoring error is the accumulation of the individual image classification errors [44][45][46].In addition, the time-series-based method avoids manipulation of a single image, thereby reducing the number of intermediate links with possible errors.Second, the time series similarity measurement can be used to identify impervious surface expansion finely, but it relies heavily on the change samples, because it is necessary to provide change samples in different time periods.The proposed approach in this study only requires samples with a certain change feature to determine the type of change at any time point.Third, the method of trajectory classification can indeed be employed to detect impervious surface expansion at different time points based on a small number of samples, but the accuracy is low due to noise because the original time series are used.In this study, noise caused by, e.g., clouds, shadows, and interannual difference, was eliminated, and only the trend trajectory of time series was preserved by segmentation fitting, thereby increasing the impervious surface expansion monitoring accuracy.

Time Series Indices Selection
Impervious surface expansion monitoring is mainly based on the change features of image time series.Therefore, the selection of indices is very important for constructing time series, and it is necessary to find indices that can highlight impervious surface expansion.In Landsat image time series change detection, indices such as the NDVI, normalized difference built-up index (NDBI), normalized burn ratio (NBR), Tasseled Cap brightness (TCB), Tasseled Cap greenness (TCG), and Tasseled Cap wetness (TCW) are often used (Table 4) [47][48][49].Using the same data set, we ran the six abovementioned indices sequentially in LandTrendr to see the trend features of the impervious surface expansion sample points corresponding to different indices and then selected the indices with relatively strong sensitivity to impervious surface expansion.In Figure 11, it can be seen that for the TCB, TCG, TCW, and NDBI, the time series trajectories of the eight sample points have different shapes and do not show any obvious trends, indicating that these indices are less sensitive to impervious surface expansion.The time series trajectories of the NDVI and NBR clearly exhibit trends and accurately monitored the change time, showing strong sensitivity to impervious surface expansion.In addition, the NDVI can effectively distinguish surface developed areas from vegetation coverage areas [49], and it has a significant effect on the detection of land use changes caused by human events.Therefore, the NDVI was used as a time series index for impervious surface expansion in this study.
these indices are less sensitive to impervious surface expansion.The time series trajectories of the NDVI and NBR clearly exhibit trends and accurately monitored the change time, showing strong sensitivity to impervious surface expansion.In addition, the NDVI can effectively distinguish surface developed areas from vegetation coverage areas [49], and it has a significant effect on the detection of land use changes caused by human events.Therefore, the NDVI was used as a time series index for impervious surface expansion in this study.

Factors Affecting the Accuracy of Impervious Surface Expansion Time
In this study, the expansion time was identified by detecting the largest descending segment of the time series.After segmentation and fitting, the redundant vertices of the original time series were removed, leaving only the vertices representing the trajectory.The impervious surface expansion time was defined as the year following the beginning of the NDVI decrease in the segment in which the time series trajectory had the largest decline.Some factors may affect this method when monitoring impervious surface expansion time.
It can be seen from the confusion matrix (Table 3) that the user accuracy was less than 90% in 2006, 2007, and 2017 but greater than 90% in the other years.A possible reason for this result is that the method has low accuracy for the impervious surface expansion time at the beginning and end of the time series and relatively high accuracy in the intermediate part of the time series.This effect can be mitigated by increasing the length of the time series.
The distribution of samples is shown in Figure 12. Figure 13a,b depict the time series of samples with the time correctly and incorrectly monitored, respectively.Comparison of these images demonstrates that the time series of the correct samples mostly had three segments.Only the time series in which impervious surface expansion occurred at the beginning and end had two segments, because there was no reference data before and after these time points.The time series of the wrong samples mostly had one or two segments.It can be considered to some extent that the shape of the time series after segmentation affected the time accuracy.When three segments were presented after segmentation, the time accuracy of the impervious surface was high; in other cases, the accuracy was low.In Figure 13, the samples with time monitoring errors show longer periods of impervious surface expansion.Five samples with time monitoring errors in 2017 exhibit short expansion periods because they were at the end of the time series.Among the remaining 13 incorrectly segmented samples, nine display expansion periods of four years or more.Thus, the impervious surface expansion period had a certain effect on the time accuracy.The longer the expansion period, the lower the time accuracy.
Mixed pixels are factors influencing precision that cannot be ignored.Because Landsat data were used in this study, the spatial resolution was 30 m, which inevitably produces some mixed pixels.For the samples with incorrect monitoring of the impervious surface expansion time in Figure 12, the Landsat and Google Earth high-resolution images of the samples were compared one-by-one, and it was found that there were mixed pixels in the misidentified samples.We selected one representative sample, as shown in Figure 14.As can be seen from the high-resolution image in Figure 14c, the actual expansion time of the sample was 2011.The expansion time monitored in this study was 2012, which is shown in Figure 14a.It can be seen that if the pixel has an impervious surface mixed with vegetation, the monitored change time may be delayed compared with the actual change time due to the influence of vegetation.Mixed pixels are factors influencing precision that cannot be ignored.Because Landsat data were used in this study, the spatial resolution was 30 m, which inevitably produces some mixed pixels.For the samples with incorrect monitoring of the impervious surface expansion time in Figure 12, the Landsat and Google Earth high-resolution images of the samples were compared one-by-one, and it was found that there were mixed pixels in the misidentified samples.We selected one representative sample, as shown in Figure 14.As can be seen from the high-resolution image in Figure 14c, the actual expansion time of the sample was 2011.The expansion time monitored in this study was 2012, which is shown in Figure 14a.It can be seen that if the pixel has an impervious surface mixed with vegetation, the monitored change time may be delayed compared with the actual change time due to the influence of vegetation.

Conclusions
Based on the NDVI time series images, LandTrendr was employed in this study to remove noise from time series, which were then segmented, fitted, and further reconstructed to obtain the spectral change trend of each pixel.Based on the analysis of the trajectory features of the NDVI time series, the trajectory features were extracted and the time of impervious surface expansion was calculated.We used a multi-level classification method to identify the scope of impervious surface expansion, thereby realizing dynamic monitoring of impervious surface expansion.The spatiotemporal evolution features of impervious surface expansion in Xinbei District, Changzhou were further analyzed.The main conclusions of this study are as follows.
(1) The time series segmentation provided accurate results in impervious surface expansion trajectory feature extraction.The LandTrendr algorithm was used to denoise, segment, fit, and further reconstruct the NDVI time series, which simplified the time series trajectory considerably

Conclusions
Based on the NDVI time series images, LandTrendr was employed in this study to remove noise from time series, which were then segmented, fitted, and further reconstructed to obtain the spectral change trend of each pixel.Based on the analysis of the trajectory features of the NDVI time series, the trajectory features were extracted and the time of impervious surface expansion was calculated.We used a multi-level classification method to identify the scope of impervious surface expansion, thereby realizing dynamic monitoring of impervious surface expansion.The spatiotemporal evolution features of impervious surface expansion in Xinbei District, Changzhou were further analyzed.The main conclusions of this study are as follows.
(1) The time series segmentation provided accurate results in impervious surface expansion trajectory feature extraction.The LandTrendr algorithm was used to denoise, segment, fit, and further reconstruct the NDVI time series, which simplified the time series trajectory considerably while retaining the trend information as much as possible.Finally, the trajectory features of the impervious surface expansion time series were described by three phases, four types, and nine parameters.(2) The multi-level classification method combining SVM and decision tree classification effectively identified the impervious surface expansion trajectories.We firstly used SVM classification to extract the pixels with impervious surface expansion trajectories.Since many pixels had similarities to impervious surface expansion time series trajectories but experienced no impervious surface expansion, we added a decision tree classification method and established classification rules to identify the pixels corresponding to transformation as non-impervious surfaces.The proposed method yielded an overall accuracy of 90.58% and a Kappa coefficient of 0.90.(3) The impervious surface expansion in Xinbei District varied considerably by region, and the expansion area in the past years showed a downward trend with alternating peaks and valleys.The southern part of Xinbei District, bordered by the downtown area of Changzhou City, witnessed rapid impervious surface expansion.The northern area was far from downtown Changzhou; therefore, its impervious surface expansion was not as fast as that in the southern area.After 2005, the impervious surface area in Xinbei District expanded by 3427.79 hectares, with an average annual expansion of about 270 hectares and declined overall.Since the Jiangdu-Yixing Expressway and Gangcheng Avenue were built in 2009 and 2014, respectively, the impervious surface expansion area exhibited peaks in these two years, resulting in alternating peaks and valleys.
In this study, vegetation-covered land changing into impervious surface was considered, but the specific type of vegetation-covered land, such as cultivated land, forest land, garden, or grassland, was not addressed.In addition, the impervious surface was not further divided by type.In future research, efforts should be made to determine the changes in various types of land during the expansion of specific types of impervious surfaces.

Figure 1 .
Figure 1.The location of Xinbei District.

Figure 3 .
Figure 3.The procedures for noise smoothing and segmentation in LandTrendr.

Figure 3 .
Figure 3.The procedures for noise smoothing and segmentation in LandTrendr.

Figure 4 .
Figure 4. Features of impervious surface expansion trajectories.(a-d) represent four types of impervious surface expansion trajectories.I, II, and III represent three phases of impervious surface expansion.NDVI, normalized difference vegetation index.

Figure 4 .
Figure 4. Features of impervious surface expansion trajectories.(a-d) represent four types of impervious surface expansion trajectories.I, II, and III represent three phases of impervious surface expansion.NDVI, normalized difference vegetation index.

Figure 5 .
Figure 5. Multi-level classification for time series trajectories.

Figure 5 .
Figure 5. Multi-level classification for time series trajectories.

Figure 6 .
Figure 6.Sample types and distribution for SVM classification.

Figure 6 .
Figure 6.Sample types and distribution for SVM classification.

23 Figure 7 .
Figure 7. Decision tree for screening the impervious surface expansion pixels.

Figure 7 .
Figure 7. Decision tree for screening the impervious surface expansion pixels.

Figure 8 .
Figure 8.Samples statistics.(a) Typical images of conversion from vegetation to impervious surface and vegetation to pond water.(b) Separability statistics of three parameters: NDVI upon completion of expansion, average NDVI after expansion, and NDVI upon completion of recovery.

Figure 8 .
Figure 8.Samples statistics.(a) Typical images of conversion from vegetation to impervious surface and vegetation to pond water.(b) Separability statistics of three parameters: NDVI upon completion of expansion, average NDVI after expansion, and NDVI upon completion of recovery.

Figure 9 .
Figure 9. Impervious surface expansion extraction results.(a) SVM classification results.(b) Multilevel classification results.(c) Comparison of areas I, II, III, and IV using high-resolution imagery.

Figure 9 .
Figure 9. Impervious surface expansion extraction results.(a) SVM classification results.(b) Multi-level classification results.(c) Comparison of areas I, II, III, and IV using high-resolution imagery.

23 Figure 10 .
Figure 10.Spatiotemporal distribution of impervious surface expansion from 2005 to 2017.(a) Spatial distribution of impervious surface expansion.(b) Area statistics of impervious surface expansion every year.

Figure 10 .
Figure 10.Spatiotemporal distribution of impervious surface expansion from 2005 to 2017.(a) Spatial distribution of impervious surface expansion.(b) Area statistics of impervious surface expansion every year.

Figure 11 .
Figure 11.Comparison of time series indices.The eight colors represent the eight sample points experiencing impervious surface expansion in 2013-2014, and the dotted lines indicate the detected impervious surface expansion time.

Figure 12 .
Figure 12.Spatial distribution of samples in Figure 13.

Figure 12 .
Figure 12.Spatial distribution of samples in Figure 13.

Figure 13 .
Figure 13.Trajectories of partial verification samples.Trajectories of samples with the time (a) correctly and (b) incorrectly monitored and trajectories of samples for which the monitoring result was incorrect in (c) 2006, (d) 2017, and (e) other years.

Figure 13 .
Figure 13.Trajectories of partial verification samples.Trajectories of samples with the time (a) correctly and (b) incorrectly monitored and trajectories of samples for which the monitoring result was incorrect in (c) 2006, (d) 2017, and (e) other years.

Figure 14 .
Figure 14.Effects of mixed pixels on time monitoring accuracy.(a) Sample trajectory, where the red dotted line indicates the detected change time.(b) Landsat images of the sample.(c) High-resolution images of the sample from Google Earth.

Figure 14 .
Figure 14.Effects of mixed pixels on time monitoring accuracy.(a) Sample trajectory, where the red dotted line indicates the detected change time.(b) Landsat images of the sample.(c) High-resolution images of the sample from Google Earth.

Author Contributions:
All authors contributed to this work.B.W. and Z.C. had the original idea of this study and designed the research.B.W. drafted the manuscript, which was revised by Z.C. and A.-X.Z.Y.H. and C.X. contributed to data processing and experiments.

Table 1 .
A list of remote sensing image data sources.
USGS, United States Geological Survey.

Table 1 .
A list of remote sensing image data sources.

Table 2 .
Parameters representing the time series features of impervious surface expansion.

Table 2 .
Parameters representing the time series features of impervious surface expansion.

Table 4 .
Calculation of time series indices.