Landsat-Based Land Cover Change in the Beijing-Tianjin-Tangshan Urban Agglomeration in 1990 , 2000 and 2010

Rapid urbanization dramatically changes the local environment. A hybrid classification method is designed and applied to multi-temporal Landsat images and ancillary data to obtain land cover change datasets. A support vector machine (SVM) classifier is used to classify multi-temporal Landsat Enhanced Thematic Mapper Plus (ETM+) images that were collected in 2000 at the pixel level. These images are also segmented with the mean shift method. The impervious surface is refined based on a combination of the segmented objects and the SVM classification results. The changed areas in 1990 and 2010 are determined by comparing the Thematic Mapper (TM) and ETM+ images via the re-weighted multivariate alteration detection transformation method. The TM images that were masked as changed areas in 1990 and 2000 are input into the SVM classifier. Land cover maps for 1990 and 2010 are produced by combining the unchanged area in 2000 with the new classes of the changed areas in 1990 and 2010. Land cover change has continuously accelerated since 1990. Remarkably, arable land decreased, while the impervious surface area significantly increased.


Introduction
Land cover and land use provide information that improves our understanding of the interactions between humans and the environment [1].Land cover influences the energy balance and the carbon and hydrological cycles and, thus, plays an important role in global change research [2][3][4].Land cover directly affects the physical characteristics of the land surface, such as soil moisture, albedo, temperature and transpiration, so many scientific studies require information regarding the spatial distribution and dynamic changes of land cover [5][6][7].Many projects, such as the International Geosphere Biosphere Programme (IGBP), the Global Observations of Forest Cover and Global Observations of Land Dynamics (GOFC-GOLD) and the Global Rain Forest Mapping (GRFM) project, have been proposed to understand the land cover and land cover changes [8][9][10].
During the past two decades, China's economic development has resulted in the rapid expansion of urban areas [11][12][13], and some urban agglomerations have gradually formed, such as the Beijing-Tianjin-Tangshan (BTT), Wuhan, Yangtze River Delta and Pearl River Delta urban agglomerations [14].The dramatic land cover changes around urban agglomerations cause a variety of problems, such as regional climate change, the urban heat island effect and increased greenhouse gas emissions [15][16][17][18][19].Many climate models, such as the Regional Atmospheric Modeling System (RAMS), the Penn State University/National Center for Atmospheric Research (PSU/NCAR) mesoscale model (MM5) and the Weather Research and Forecasting Model (WRF), require land cover as input data to analyze the effects of land cover changes on the climate [20][21][22][23][24]. BTT is the largest urban agglomeration in North China [14].Using climate models to simulate the effects of land cover changes on the regional climate in the BTT region requires long-term land cover products.
Currently, many land cover and land use datasets that cover the BTT region are freely available.Datasets with resolutions from 300 m-1000 m were primarily developed based on Medium Resolution Imaging Spectrometer (MERIS), Advanced Very High Resolution Radiometer (AVHRR), Moderate Resolution Imaging Spectroradiometer (MODIS) and Satellite Pour l'Observation de la Terre (SPOT) Vegetation [25][26][27][28].These datasets represent the land surface characteristics in different years, but none of these datasets are continuous from 1990-2010.Furthermore, datasets with coarse resolution have low classification accuracy, especially in urban areas [21,29,30].The information from classification products with 30-m spatial resolution is more detailed than that from coarser products, but the number of classes in classification products with 30-m resolution is generally lower than what is required by climate models [1,31,32].In addition to the above global-and national-scale datasets, many land cover classifications of Beijing have been conducted, but the other regions in the BTT urban agglomeration have rarely been addressed [33][34][35][36].
Differences between land cover products are more susceptible to classification methods and classification systems [37].The land covers must have identical semantics or extremely high accuracy to monitor any changes [38].The above-mentioned land covers have obvious differences in thematic classes, classification methods and data.Comparing inconsistent maps creates unreliable conclusions [38].Consistent land covers in 1990, 2000 and 2010 must be developed to understand the land cover changes from the expansion of the BTT urban agglomeration.
Images from remote sensing satellites have become established as the main basis for the development of land covers [26,28,39].Classifications can be divided into pixel-based and object-based classifications based on pixels or pixel clusters [40].The 'salt and pepper' problem often occurs in pixel-based classifications, so object-based classifications were developed to reduce the noise [35, [41][42][43].Segmentation is a critical step in object-based classifications.The segmentation of an image combines a group of spectrally-similar and spatially-adjacent pixels as an object [44][45][46].Generally, some parameters must be set for image segmentation, such as the spatial scale, bandwidth and shapes; however, selecting one parameter that is suitable for all land cover classes is difficult [35].Meanwhile, integrating expert knowledge and the effective attributes of the objects to be classified is time consuming [47].Therefore, methods to combine pixel-based and object-based approaches have been used in some classifications [31,[48][49][50].
Change detection is an important component of land cover studies.As a post-classification method, the comparison of classification results can be performed to detect changes [51].Furthermore, pre-classification algorithms, such as image differencing, change vector analysis, image regression and image ratios, have been used to distinguish changed areas between images that are acquired on different dates [52][53][54].Post-classification change detection could reduce the influence of radiometric calibration and sensor differences between different images, but the detection accuracy greatly depends on the classification accuracy of each land cover [55,56].In contrast, pre-classification approaches can avoid the cumulative errors of the two classification results because only the changed area is updated.Generally, pre-classification approaches require images that are acquired on similar dates, but the iteratively re-weighted multivariate alteration detection transformation (IR-MAD) algorithm can be used to combine the relative radiometric matches and change detection [57].The IR-MAD algorithm has no strict limit on acquiring images in the same phenological period.In this study, the IR-MAD algorithm is used to identify changed areas.
The aim of this study is to identify consistent land covers in 1990, 2000 and 2010 in the BTT urban agglomeration by using Landsat images.The classification system and accuracies of the land covers should satisfy the requirements of the RAMS climate model.To achieve this goal, we adopted a hybrid classification approach that integrated the advantages of both pixel-based and object-based approaches to classify multi-temporal Landsat images in 2000.The land covers in 1990 and 2010 were produced based on the land cover in 2000 by using the pre-classification change detection method.The characteristics of the urbanization in the BTT region were analyzed based on these three land covers.

Study Area
The BTT urban agglomeration (38

Data and Pre-Processing
The Landsat satellite series has continually observed Earth since 1972 and has accumulated an enormous number of time series images [38].Landsat images have been widely used in land cover classification because of their stable imaging quality [1,31,32,[58][59][60].In this study, Level 1 terrain-corrected (L1T) Landsat images were selected as the primary data source for the land cover classifications.The cloud coverage in the selected Landsat images was less than 5%.The Landsat images were atmospherically corrected to produce the surface reflectance by using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS).LEDAPS uses the 6SV (Second Simulation of a Satellite Signal in the Solar Spectrum vector code) radiative transfer code.The aerosol characterizations are derived independently from each Landsat acquisition, assume a fixed continental aerosol type and use ancillary water vapor [61][62][63].The L1T Landsat images were processed to surface reflectance products by applying LEDAPS.The influence of the difference in the illumination geometry in the Thematic Mapper (TM)/ETM+ images was substantially eliminated.
Covering the entire study area required TM/ETM+ images from six scenes.The images that were used in this study are listed in Table 1.The images that were acquired in 1990 and 2010 were Landsat 5 TM images, and those that were acquired in 2000 were Landsat 7 ETM+ images.The data quality of Landsat 7 ETM+ images is superior to that of Landsat 5 TM images, so three Landsat 7

Data and Pre-Processing
The Landsat satellite series has continually observed Earth since 1972 and has accumulated an enormous number of time series images [38].Landsat images have been widely used in land cover classification because of their stable imaging quality [1,31,32,[58][59][60].In this study, Level 1 terrain-corrected (L1T) Landsat images were selected as the primary data source for the land cover classifications.The cloud coverage in the selected Landsat images was less than 5%.The Landsat images were atmospherically corrected to produce the surface reflectance by using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS).LEDAPS uses the 6SV (Second Simulation of a Satellite Signal in the Solar Spectrum vector code) radiative transfer code.The aerosol characterizations are derived independently from each Landsat acquisition, assume a fixed continental aerosol type and use ancillary water vapor [61][62][63].The L1T Landsat images were processed to surface reflectance products by applying LEDAPS.The influence of the difference in the illumination geometry in the Thematic Mapper (TM)/ETM+ images was substantially eliminated.
Covering the entire study area required TM/ETM+ images from six scenes.The images that were used in this study are listed in Table 1.The images that were acquired in 1990 and 2010 were Landsat 5 TM images, and those that were acquired in 2000 were Landsat 7 ETM+ images.The data quality of Landsat 7 ETM+ images is superior to that of Landsat 5 TM images, so three Landsat 7 ETM+ images were used for every scene in 2000 [60,64].Multi-temporal ETM+ images were utilized in the classification for 2000 to reduce spectral confusion.The scene selection was based on the vegetation dynamics of the target land cover types over a growing season.Crops in the study area could be harvested twice in one year.However, some arable lands were cultivated only once in one year.Multi-temporal images can compensate for cropland classification errors from harvesting or an absence of planting and capture agricultural changes during the crop's life [65].Multi-temporal images also make full use of phenological differences in different classes.The three acquisition dates of the ETM+ images covered the early, peak and late growing seasons [66].The choice of appropriate acquisition dates would facilitate maximizing the differences between types.The major green vegetation in the BTT region was agricultural crop.Considering the phenology of the crops and the maximum separability among dry farmland, paddy field and forest, the three acquired dates were in March, August and October.We experienced difficulty selecting three temporal images in one year because of the influences of clouds and rain.Therefore, high-quality images that were acquired in the previous or following year were selected.One Landsat 5 TM image was selected for every scene in 1990 and 2010.An extra TM image was added for scene 121/33 in 1990 to compensate for the influence of cloud coverage.In addition to the Landsat images, which constituted the primary data source, Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) data were included in the classification [67].Forest, shrub and grass were easily confused with agricultural crops in the study area because of the spectral similarity in the Landsat images.However, these forest, shrub and grass were mainly distributed in mountainous areas, and agricultural crops were mainly distributed in plains in the study area.A large difference in slope was observed between the two terrains.In addition, the Conversion of Cropland to Forest Program (CCFP) has been implemented by the Chinese government since 1999 [68].In the CCFP, cultivated land along steep slopes is converted into forest, shrub and grass [69].Therefore, any forest, shrub and grass on mountains exhibits obvious slope differences between cultivated lands.Raw three-arc-second-resolution DEM data were mosaicked and re-projected to the same projection as the TM/ETM+ images and then cut to the size same as each scene.The DEM was resampled to 30-m horizontal resolution by using the nearest neighbor resampling method.The Landsat images and the slopes from the DEM were used as the input data for the SVM classifier.Only the Landsat images were used when using the SVM to train the samples, and the accuracy of 5-fold cross-validation was 87.98%.When the slope data were added, the accuracy reached 92.63%.

Methods
The study area included many rural settlements that belonged to the impervious surface class.During the growing season, these rural settlements were extensively sheltered by trees that were planted in streets or courtyards.During the dormant season, the spectral characteristics of the rural settlements were similar to those of the bare arable land.This complicated environment reduced the classification accuracy of rural settlements when using pixel-based classification methods, and rural settlements were typically subject to omission errors.The segment object was the cluster of pixels of various classes.The characteristics that described a single pixel were not suitable to the object.Using the SVM to classify the segment objects required more characteristics.Image segmentation was performed and combined with the SVM classification results to improve the classification accuracy of the impervious surface class [31,70].
The flow diagram of the classification process is shown in Figure 2.During classification, the multi-temporal ETM+ images from 2000 and the slope data were first classified scene-by-scene by the SVM classifier.The pixel-based and object-based approaches were integrated to improve the accuracy of the impervious surface.Therefore, the mean shift algorithm was used to segment the three multi-temporal Landsat ETM+ images to generate objects.The segmented objects were combined with the SVM classification result to revise the impervious surface class.The pre-classification change detection method was used to find the changed area; thus, one of the three images from 2000 was selected for comparison with images from 1990 and 2010 by using the IR-MAD approach [52,57].Only the areas on the Landsat images from 1990 and 2010 that were masked by the changed area were classified with the SVM classifier.The land cover for 2000 was updated with the classification results of the changed area to develop the land cover for 1990 and 2010.

Classification System and Training Samples
Considering the spatial resolution of the TM/ETM+ images and the requirements of some climate models, the classification system comprised 9 classes that referenced some existing classification systems: forest, shrub, grass, dry farmland, paddy field, bare land, impervious surface, aquatic vegetation and water bodies.
The training samples were primarily collected by visual interpretation of the ETM+ images from 2000 with reference to Google Earth.The rest were collected from field observations in 2000.The sample points that were collected by visual interpretation were first randomly distributed in the six scenes.Then, a polygon that included pixels with the same class was manually drawn at a single point location.If the random sample was not in a pure pixel cluster, we selected a nearby polygon with the same class as the random points.Some land cover classes, such as water bodies, bare land and paddy fields, had lower area proportions in the study area.If only random samples were used, the sample counts of these classes would be significantly lower than those of the other classes.More samples were added by artificial selection for the classes with lower area proportions to mitigate the influence of these unbalanced sample counts.Finally, 327 polygons (26893 pixels) were identified.

Classification System and Training Samples
Considering the spatial resolution of the TM/ETM+ images and the requirements of some climate models, the classification system comprised 9 classes that referenced some existing classification systems: forest, shrub, grass, dry farmland, paddy field, bare land, impervious surface, aquatic vegetation and water bodies.
The training samples were primarily collected by visual interpretation of the ETM+ images from 2000 with reference to Google Earth.The rest were collected from field observations in 2000.The sample points that were collected by visual interpretation were first randomly distributed in the six scenes.Then, a polygon that included pixels with the same class was manually drawn at a single point location.If the random sample was not in a pure pixel cluster, we selected a nearby polygon with the same class as the random points.Some land cover classes, such as water bodies, bare land and paddy fields, had lower area proportions in the study area.If only random samples were used, the sample counts of these classes would be significantly lower than those of the other classes.More samples were added by artificial selection for the classes with lower area proportions to mitigate the influence of these unbalanced sample counts.Finally, 327 polygons (26893 pixels) were identified.
The samples were identified by visual interpretation and drawing polygons, so the classes of samples may have included mistakes.All of the samples were input into the SVM classifier that was trained with the same samples to check for mistakes.If the predicted class was not the same as that input, the class of the sample was subjected to additional checks.If the class of the sample was correct according to Google Earth images, the class was retained, even if the SVM classifier output a different class.If the class of the sample was incorrect, the class was revised.This process was repeated until we were confident that the samples included no mistakes.

SVM Classification
Recently, the SVM classifier has been widely used in remote sensing classification with good results [71][72][73].SVM uses the kernel function to map the spectral feature space to high-dimensional space and calculates the inner product in this high-dimensional space.SVM performs classification by using an optimal hyperplane, which is determined based on the limited supporting samples [74].SVM can achieve high classification accuracy with relatively small numbers of training samples [73].In this study, the Libsvm (Version 3.18) software package was used [75].This software package can normalize the input data and search for the optimal parameters.Bands 1-5 and 7 of the Landsat images were set as the main input data for the SVM classifier.The slope data from the DEM were also included in the classification.Therefore, the input data for the year 2000 had 19 layers, including one slope layer and 18 bands (each scene corresponded to 6 bands) from the Landsat images of three scenes.Each scene was individually classified, and the results of six scenes were mosaicked together.
During classification, the 19 layers were normalized by using the scale function, and the radial basis function (RBF) was set as the kernel function.The Libsvm software package performed numerous classification tests to search for the optimal parameters by setting the searching step length of the parameters.Specifically, the samples of every scene were first scaled to the range [−1, +1].Then, svm_type was set as c_svc; the 10-fold cross-validation was performed; and the step lengths of the penalty parameter and the gamma parameter in RBF were set to 1.The iteration termination threshold was 0.001.SVM classification was conducted once the optimal model parameters were determined.

Image Segmentation
The mean shift algorithm has been widely used in image segmentation and identification tracking [76,77].As a kernel density estimate method, the mean shift does not require prior knowledge of the number of segmentations and sets no limits on the segmented shape [78].In this study, the EDISON V1.1 software was used.The 18 reflective bands of the three scenes of Landsat images were transformed based on principal component analysis (PCA) to reduce the calculation iterations of the mean shift.The first three principal components were selected as the segmented layers.The Gaussian function was taken as the kernel function, and the bandwidths of the geometry space and color space were set to 15 and 9.5, respectively.
In this study, the segmented object was not set as the input data for the SVM classifier.A set of appropriate thresholds was used according to the area percentages of the classes within an object to determine whether the object belonged to the impervious surface class.Forest and dry farmland were the main factors because the interference factors to the rural settlements mainly included the trees and dry farmland around or inside the rural settlements.The segmented objects were overlaid with the SVM classification to refine the impervious surface.Generally, if the area fractions of impervious surface (A i ), forest (A f ) and dry farmland (A d ) within an object met the judgment conditions, the polygon was designated as an impervious surface.These thresholds were selected based on multiple tests.The following rules were designed through the stratified extraction of the impervious surface: If 40% ≤ A i < 50% and 40% ≤ A f < 50%, class = impervious surface, If 30% ≤ A i < 40% and 30% ≤ A f < 40% and 10% ≤ A d < 30%, class = impervious surface (3) where A i , A f and Ad are the area fractions of impervious surface, forest and dry farmland, respectively.The rural settlements were mainly subjected to omission errors, so objects were classified as impervious surface when the Ai pixels were greater than 50% in a segmented object.If an object had more than 40% impervious surface pixels and forest pixels, the total pixels of these two classes were more than 80%.Basically, the rural settlements matched this situation.The majority of impervious surface objects were identified through these two equations.The proposed Equation (3) considered that the proportions of rural settlements, trees and farmland were all relatively low.
Furthermore, the area percentage of the small impervious surface patch may have been less than 30% if a small impervious surface patch was located in a large object.Therefore, the impervious surface patch would be deleted, even if the impervious surface patch was correctly classified.An extra rule was set to avoid this situation: if one impervious surface patch comprised more than 10 adjacent pixels, the impervious surface patch would be reserved regardless of the percentages of the classes in one object.If an impervious surface patch was deleted, the adjacent class with the largest area would be assigned to the deleted location.

Land cover Change
The IR-MAD algorithm was used to detect the changed area, and then, the TM images of the changed area were classified by the SVM classifier.IR-MAD judged the unchanged pixels by iteratively calculating the chi-square distribution of the difference in the bi-temporal images [52,57].The IR-MAD linearly transformed the original images R and T to new images U and V based on the canonical correlation analysis (CCA) in Equation ( 4): where U and V are the canonical correlation images, a T and b T are the transform matrixes and R and T are the bi-temporal images.
The bands of the canonical correlation images were referred to as the canonical variates.The canonical variates were arranged according to the correlations.If the correlation was relatively low, this pair of canonical variates had more change information.The differences in each pair of canonical variables were mutually uncorrelated.
According to the central limit theorem, the differences in the canonical variates approximately corresponded to a Gaussian distribution, and the sum of squares of the differences in the canonical variates corresponded to a chi-square distribution [52,57].The differences in the canonical variates and the weights of unchanged pixels were calculated by using Equations ( 5)-( 7): Pr (no change) = 1 − P χ 2 ;N (Z), (7) where i is the band of the image, M i is the difference in the variates, σ M i is the standard deviation of M i , N is the total number of bands, P χ 2 ;N (Z) is the χ 2 distributed with N degrees of freedom and Pr (no change) is the weight.
Pr (no change) was used to weigh each pixel.Then, Equations ( 4)-( 7) were iterated until no significant changes were observed in the canonical correlations.These pixels could then be determined as unchanged pixels if their chi-square distribution probability was lower than 0.9 [52].
In this study, we calculated the chi-square images from 2010 and 1990 based on the images from 2000.The threshold was slightly lower than 0.95, which was set by the developer, so some unchanged pixels were selected in the identified area.Although some extra unchanged pixels were selected, mostly changed pixels were chosen, so the omission error was reduced.
The samples were collected from the unchanged area to classify the changed pixels with the SVM classifier.The land covers from 1990 and 2010 were produced by updating the land cover from 2000 with the classification of the changed areas for 1990 and 2010.Figure 3 shows this process.

Accuracy Assessment
The assessment samples were selected through stratified random sampling.Bare land and aquatic vegetation had the smallest area proportions, so we supplemented some points to these two classes.The class of the assessment sample was identified by visual interpretation.A total of 526 points from these nine classes were collected to assess the classification accuracy.The error matrix was constructed to calculate the overall accuracy, Kappa coefficient, user's accuracy and producer's accuracy for the three periods [79].The overall accuracy is the proportion of the area that was correctly mapped.The user's accuracy is the proportion of correctly-classified pixels with regard to all of the pixels that are classified as this class in the classified image.The producer's accuracy is the proportion of correctly-classified pixels with regard to all of the pixels of that ground truth class [80].The overall accuracy, user's accuracy and producer's accuracy were adjusted in terms of the unbiased estimator of the area proportions [80].The specific calculation process is shown in Equations ( 8)-( 11) [81].The results of the accuracy assessment are shown in Table 2.

  
where .i N is the number of pixels of class i in the entire map, N is the number of pixels in the entire map, .i n is the number of pixels from a sample of class i and ij n is the sample count in row i (map category) and column j (reference category) in the error matrix of the sample counts.

Accuracy Assessment
The assessment samples were selected through stratified random sampling.Bare land and aquatic vegetation had the smallest area proportions, so we supplemented some points to these two classes.The class of the assessment sample was identified by visual interpretation.A total of 526 points from these nine classes were collected to assess the classification accuracy.The raw assessment points and error matrixes were shown in supplementary materials.The error matrix was used to calculate the overall accuracy, Kappa coefficient, user's accuracy and producer's accuracy for the three periods [79].The overall accuracy is the proportion of the area that was correctly mapped.The user's accuracy is the proportion of correctly-classified pixels with regard to all of the pixels that are classified as this class in the classified image.The producer's accuracy is the proportion of correctly-classified pixels with regard to all of the pixels of that ground truth class [80].The overall accuracy, user's accuracy and producer's accuracy were adjusted in terms of the unbiased estimator of the area proportions [80].The specific calculation process is shown in Equations ( 8)-( 11) [81].The results of the accuracy assessment are shown in Table 2 where N i. is the number of pixels of class i in the entire map, N is the number of pixels in the entire map, n i. is the number of pixels from a sample of class i and n ij is the sample count in row i (map category) and column j (reference category) in the error matrix of the sample counts.pPA = pij p .j (9) pUA = pij ISPRS Int.J. Geo-Inf.2017, 6, 59 9 of 21 pOA = ∑ pkk (11) where pPA is the corrected producer's accuracy, pUA is the corrected user's accuracy and pOA is the corrected overall accuracy.p .j is the sum of the pij in column j of the estimated error matrix.p i. is the sum of the pij in row i of the estimated error matrix.pkk is the diagonal value of the estimated error matrix.The land cover in 2000 had the highest overall accuracy and Kappa coefficient, which were 89.13% and 0.87 respectively.The overall accuracies and Kappa coefficients were 87.94% and 87.76% and 0.85 and 0.85 for 1990 and 2010, respectively.Overall, the accuracy for 2000 was higher than that for 1990, and the accuracy for 1990 was higher than that for 2010; however, the difference was minor.We were more concerned about the dry farmland, forest and impervious surface in the nine classes.Indeed, these three classes comprised approximately 70% of the study area and were more seriously affected by urbanization.The user's accuracy and producer's accuracy of these three classes all exceeded 80%.
A sub-area was selected to compare the land cover in 2000 with China's Land Use/cover Dataset (CLUDs) for the year 2000 [32].This sub-area includes mountains and plains in Landsat Path/Row 123/032.The land covers that were classified in this study and the CLUDs are shown in Figure 4. Figure 4a-c shows the Landsat images in 1990, 2000 and 2010.Figure 4d-f shows the land covers that were classified in this study for Figure 4a-c.Figure 4g is the CLUD for the year 2000.This CLUD was derived by visual interpretation [32].When comparing the land cover in 2000 with the CLUD for the year 2000, the main differences were the distribution of grass and shrub.The interpreted boundaries of grass and shrub in this CLUD were artificial drawings.These polygons, which were drawn by humans, could not adequately present the interleaving and transitions among vegetation.However, if no water was present in a riverbed, this area would be classified as bare land or grass in the land cover for the year 2000.In the CLUD for the year 2000, this riverbed was interpreted as a river.This difference shows the flexibility in visual interpretation.

Impervious Surface Refinement
The impervious surface class from the SVM included many incorrect trivial speckles because this class is affected by the complexity of the environment.However, integration with segmented objects effectively reduced these incorrect trivial speckles.
The impervious surface of a sub-image was visually interpreted with reference to the high-resolution images to evaluate the effect of the impervious surface refinement.Detailed information is shown in Table 3.In this table, 'no change' means that no changes occurred; 'add correct' indicates that the impervious surface was correctly added to the SVM classification; 'add error' indicates that the impervious surface was incorrectly added; and 'remove correct' and 'remove error' refer to correctly or incorrectly deleted impervious surface, respectively.The percentage was calculated based on the area of impervious land that was identified by the SVM.Overall, the percentage of correct refinement far exceeded the percentage of incorrect refinement.In this sub-area, the percentages of impervious surface pixels that were added and removed correctly were 18.99%

Impervious Surface Refinement
The impervious surface class from the SVM included many incorrect trivial speckles because this class is affected by the complexity of the environment.However, integration with segmented objects effectively reduced these incorrect trivial speckles.
The impervious surface of a sub-image was visually interpreted with reference to the high-resolution images to evaluate the effect of the impervious surface refinement.Detailed information is shown in Table 3.In this table, 'no change' means that no changes occurred; 'add correct' indicates that the impervious surface was correctly added to the SVM classification; 'add error' indicates that the impervious surface was incorrectly added; and 'remove correct' and 'remove error' refer to correctly or incorrectly deleted impervious surface, respectively.The percentage was calculated based on the area of impervious land that was identified by the SVM.Overall, the percentage of correct refinement far exceeded the percentage of incorrect refinement.In this sub-area, the percentages of impervious surface pixels that were added and removed correctly were 18.99% and 13.48%, respectively.The corresponding percentages for incorrect addition and removal were 2.89% and 1.99%, respectively.Therefore, the accuracy of the impervious surface classification was improved by 27.59%.Figure 5 shows the local impervious surface that was classified by the SVM, the objects that were segmented by the mean shift method and the refined impervious surface after synthesizing the SVM and mean shift results.Figure 5a is a false-color Landsat image.Figure 5b shows the vector boundary of the segmented objects when overlaid onto the sub-image in Figure 5a.The outlines of the impervious surface are clearly delineated.In Figure 5c, the impervious surface from the SVM contains many trivial speckles.Figure 5d shows the refined impervious surface.Detailed information regarding the integration is presented in Figure 5e.
ISPRS Int.J. Geo-Inf.2017, 6, 59 11 of 21 and 13.48%, respectively.The corresponding percentages for incorrect addition and removal were 2.89% and 1.99%, respectively.Therefore, the accuracy of the impervious surface classification was improved by 27.59%.Figure 5 shows the local impervious surface that was classified by the SVM, the objects that were segmented by the mean shift method and the refined impervious surface after synthesizing the SVM and mean shift results.Figure 5a is a false-color Landsat image.Figure 5b shows the vector boundary of the segmented objects when overlaid onto the sub-image in Figure 5a.The outlines of the impervious surface are clearly delineated.In Figure 5c, the impervious surface from the SVM contains many trivial speckles.Figure 5d shows the refined impervious surface.Detailed information regarding the integration is presented in Figure 5e.

Change Mask
A threshold of 0.9 was manually chosen for the chi-square image based on the IR-MAD algorithm to extract the changed area.Figure 6 shows the sub-images, which contain complex changes in various classes.Figure 6a,b is the false-color Landsat sub-image that were acquired on 2 August 1999 and 22 September 2009, respectively.The main classes in the sub-images were an airport, towns, rural settlements, dry farmland and rivers.The airport site at the center grew almost two-fold from 1999-2009.The bright cyan areas correspond to buildings and roads.Forests are shown in dark red in the upper corner along the river.A substantial amount of dry farmland (bright red or light cyan) became impervious surface.Figure 6c is the chi-square image; the grayscale values reflect the intensity of the changes.The chi-square image exhibited a clear response to these class changes.Figure 6d is the mask of the changed area in yellow, and the background is the sub-image in Figure 6a.Clearly, dramatic changes occurred around the airport.

Change Mask
A threshold of 0.9 was manually chosen for the chi-square image based on the IR-MAD algorithm to extract the changed area.Figure 6 shows the sub-images, which contain complex changes in various classes.Figure 6a,b is the false-color Landsat sub-image that were acquired on 2 August 1999 and 22 September 2009, respectively.The main classes in the sub-images were an airport, towns, rural settlements, dry farmland and rivers.The airport site at the center grew almost two-fold from 1999-2009.The bright cyan areas correspond to buildings and roads.Forests are shown in dark red in the upper corner along the river.A substantial amount of dry farmland (bright red or light cyan) became impervious surface.Figure 6c is the chi-square image; the grayscale values reflect the intensity of the changes.The chi-square image exhibited a clear response to these class changes.Figure 6d is the mask of the changed area in yellow, and the background is the sub-image in Figure 6a.Clearly, dramatic changes occurred around the airport.

Land Cover and Change
The land cover in the BTT region has dramatically changed over the past two decades.The total areas of the study region in 1990, 2000 and 2010 were 55,508 km 2 , 55,660 km 2 and 55,855 km 2 , respectively.Thus, the area increased by 152 km 2 and 195 km 2 from 1990-2000 and from 2000-2010, respectively.The increased areas corresponded to sea reclamation and thus were concentrated in the coastal region.The increased land from the sea mainly included impervious surface and bare land.The land cover maps for the three periods are shown in Figure 7, and the local details are shown in Figure 4.

Land Cover and Change
The land cover in the BTT region has dramatically changed over the past two decades.The total areas of the study region in 1990, 2000 and 2010 were 55,508 km 2 , 55,660 km 2 and 55,855 km 2 , respectively.Thus, the area increased by 152 km 2 and 195 km 2 from 1990-2000 and from 2000-2010, respectively.The increased areas corresponded to sea reclamation and thus were concentrated in the coastal region.The increased land from the sea mainly included impervious surface and bare land.The land cover maps for the three periods are shown in Figure 7, and the local details are shown in Figure 4.The areas and proportions of the classes were adjusted to correct for the estimated error based on the assessment samples.The adjusted areas and proportions are shown in Table 4.The proportions of each class for the three periods are also shown in Figure 8.The areas and proportions of the classes were adjusted to correct for the estimated error based on the assessment samples.The adjusted areas and proportions are shown in Table 4.The proportions of each class for the three periods are also shown in Figure 8.The sum of the areas of dry farmland, forest and impervious surface always exceeded 70% of the total study area during the two analyzed decades.The most remarkable changes were the reduction in arable land and the increase in impervious surface between 1990 and 2010.Dry farmland changed the most significantly from 1990-2000, shrinking by 1260 km 2 , with its proportion of the total The sum of the areas of dry farmland, forest and impervious surface always exceeded 70% of the total study area during the two analyzed decades.The most remarkable changes were the reduction in arable land and the increase in impervious surface between 1990 and 2010.Dry farmland changed the most significantly from 1990-2000, shrinking by 1260 km 2 , with its proportion of the total area decreasing from 47.1% down to 44.7%.Subsequently, dry farmland decreased further to 42.09% from 2000-2010.The area of paddy fields was approximately 5.76% in 1990 and decreased to 4.92% and 2.95% in 2000 and 2010, respectively.In contrast, the area of impervious surface increased the most during the studied period.The area of impervious surface increased from 8.72%-11.45%from 1990-2000 and then to 16.94% by 2010, nearly doubling in only two decades.The increased area of impervious surface was mostly attributable to the conversion of arable land.Indeed, approximately 90% of the increased area of impervious surface between 1990 and 2000 was originally arable land; from 2000-2010, this value was approximately 80%.
Compared to the changes seen in the areas of dry farmland and impervious surface, the changes in the areas of other classes were relatively small.The forest area exhibited a slight decrease from 1990-2000, followed by an increase of 0.72% by 2010.The decreased forest area was mainly converted to dry farmland and impervious surface from 1990-2000.Beginning in 1999, China began to convert cropland to forest, so some dry farmland with steep slopes and some shrub land were converted to forest land.Although water bodies had some increments (i.e., the enclosed sea), the total area that was covered by water bodies continuously decreased.Clearly, the water levels of the reservoirs in this area declined.

Integration of Pixel-Based and Object-Based Classification Methods
The combination of pixel-based classification and object-based image segmentation was critical in this study.The land surface characteristics of the study area are very complex.For example, some agricultural land was not planted, so its spectral characteristics were similar to those of bare land and buildings.Many greenhouses were located on agricultural land and appeared similar to rural settlements in the multi-spectral images.Furthermore, buildings in some rural settlements were sheltered by trees.This situation is not as prevalent in urban areas because buildings are generally higher than trees.During classification, some rural settlements were classified as several disconnected components when using the SVM classifier.However, the buildings and trees in these rural settlements could be segmented as one object.These objects could comprise one or several types of land cover pixels, so the spectral characteristics of these objects were more complex than those of single pixels.Therefore, the classification of these objects required more information, such as the shape, perimeter, area and location relative to other objects.When we attempted to use some object features, such as the mean, median, maximum and minimum, as the input data for the SVM classifier, the classification results were not improved relative to those that were only based on pixels.Indeed, the mean, maximum and minimum values of the objects widely varied, especially for objects with pixels of several classes.Thus, the identification of objects requires more knowledge, such as the relationship between objects, the distribution pattern of pixels in an object and the shape of the object.Substantial additional work is required to explore the knowledge for every class to classify objects with the SVM classifier.
Image segmentation generally requires segmentation parameters to be set.However, large differences in spatial scale existed between the land cover classes, so properly segmenting the images for all of the land cover classes with one spatial scale parameter was difficult.In this study, the segmentation parameters were set to a scale that was suitable for rural settlements to improve the classification accuracy of the impervious surface.We mainly considered the proportions of the buildings, trees and agricultural land in a single object to combine the advantages of both pixel-based and object-based approaches.In this study, the error that was associated with impervious surface based on the SVM classifier was mainly omission error.Therefore, the thresholds were set to supplement the impervious surface when combining the SVM classification results and segmented objects.Finally, the gaps in impervious surface were filled, and the small misclassified patches were removed.The integrating rules could be easily understood, but the rules were constructed by many experiments.These rules were empirical rules and should be redesigned for different regions or different types.In addition, types with large-scale differences will not be adequately segmented because of the limitations of the segmentation-scale parameters.Self-adaptive segmentation could be a potential solution.

Change Detection
Using a consistent classification method is critical to ensure that the classification results for different periods are comparable.The spectral characteristics of different land covers usually have some confusion or overlap and are not always separated by a clear-cut boundary.The transitions between some classes, such as between bare land and sparse grassland and between wetlands and water bodies, are continuous.Classification errors often appear near the boundaries between these classes because of their spectral similarity.Each satellite image has unique Sun-illumination geometry and weather characteristics.Despite using the same classifier to classify images that are acquired at different dates, ensuring that the boundaries of the unchanged land cover in the transition region are similarly identified is difficult.This type of discrepancy in boundary-adjacent areas can result in errors when comparing land covers between different periods.Generally, only a portion of the land cover in a region changes over a given period.The error in land cover changes is lower if only the changed areas are updated and the unchanged areas are retained.However, updating the changed area particularly relies on the ability to accurately identify the changed area.If the identified changed area is substantially greater than the actual changed area, the effect of updating is the same as classifying the entire image.In contrast, the omission error increases if the identified changed area is less than the actual changed area.If the changed area can be appropriately identified, however, updating only the changed area ensures that the new land cover product is consistent with the base product.
In this study, the threshold of changed pixels was set slightly lower than the recommended threshold by the software developer, so the identified changed area was slightly greater than the actual changed area.This choice was made to reduce the omission error.Overall, the method of updating the changed area reduced the required computation for the classification and improved the consistency of the land covers between different periods.

Conclusions
China's urbanization process has continuously progressed, and some urban agglomerations have formed.This urbanization has led to more frequent and intense land cover changes, thus affecting the local climate.Multi-temporal Landsat images and ancillary data were classified to extract LCLUC information for the BTT urban agglomeration.
In this study, the SVM classifier was used to classify multi-temporal Landsat images that were acquired in 2000.The objects that were identified via image segmentation were integrated with the SVM land cover classification results to improve the classification accuracy of impervious surface.The change detection method IR-MAD was used to update the land cover in 1990 and 2010 based on the land cover in 2000.Finally, consistent land covers were developed for 1990, 2000 and 2010.Remarkably, the area of arable land in the BTT urban agglomeration decreased over the past two decades, while the impervious surface area increased.The area of impervious surface increased nearly two-fold from 1990-2010.Meanwhile, the decreased area of arable land, including dry farmland and paddy fields, was slightly less than the increased area of impervious surface.The method that was used in this study was proven to be effective after validating the classification results.Thus, this method could be used for similar applications.
• 28 -41 • 3 N, 115 • 23 -119 • 48 E) includes 5 cities: Beijing, Langfang, Tianjin, Tangshan and Qinhuangdao.As the Chinese capital, Beijing promotes the economic development of its surrounding areas.Similarly, as a municipality, Tianjin represents the economic center of the Pohai Bay region.The influence of Beijing and Tianjin has served to accelerate the economic development of the adjacent cities.The location of the study area is shown in Figure 1.The study area is surrounded by the Taihang Mountains to the west, the Yanshan Mountains to the north and the Pohai Sea to the east.The main portion of the study area is arable land and is located in the northeastern area of the North China plain.Wheat and maize are widely planted in the plains areas, whereas forest and shrub are the major classes of vegetation in the mountainous areas.ISPRS Int.J. Geo-Inf.2017, The BTT urban agglomeration (38°28′-41°3′ N, 115°23′-119°48′ E) includes 5 cities: Beijing, Langfang, Tianjin, Tangshan and Qinhuangdao.As the Chinese capital, Beijing promotes the economic development of its surrounding areas.Similarly, as a municipality, Tianjin represents the economic center of the Pohai Bay region.The influence of Beijing and Tianjin has served to accelerate the economic development of the adjacent cities.The location of the study area is shown in Figure 1.The study area is surrounded by the Taihang Mountains to the west, the Yanshan Mountains to the north and the Pohai Sea to the east.The main portion of the study area is arable land and is located in the northeastern area of the North China plain.Wheat and maize are widely planted in the plains areas, whereas forest and shrub are the major classes of vegetation in the mountainous areas.

Figure 1 .
Figure 1.The map shows the location of the study area in China.The magnified map is a false-color mosaic (red-green-blue (RGB): Bands 4, 3 and 2, respectively) of Landsat Enhanced Thematic Mapper (ETM) images.The cities of Beijing, Langfang, Tianjin, Tangshan and Qinhuangdao compose the urban agglomeration.

Figure 1 .
Figure 1.The map shows the location of the study area in China.The magnified map is a false-color mosaic (red-green-blue (RGB): Bands 4, 3 and 2, respectively) of Landsat Enhanced Thematic Mapper (ETM) images.The cities of Beijing, Langfang, Tianjin, Tangshan and Qinhuangdao compose the urban agglomeration.
ISPRS Int.J. Geo-Inf.2017, 6,59 5 of 21 performed and combined with the SVM classification results to improve the classification accuracy of the impervious surface class[31,70].The flow diagram of the classification process is shown in Figure2.During classification, the multi-temporal ETM+ images from 2000 and the slope data were first classified scene-by-scene by the SVM classifier.The pixel-based and object-based approaches were integrated to improve the accuracy of the impervious surface.Therefore, the mean shift algorithm was used to segment the three multi-temporal Landsat ETM+ images to generate objects.The segmented objects were combined with the SVM classification result to revise the impervious surface class.The pre-classification change detection method was used to find the changed area; thus, one of the three images from 2000 was selected for comparison with images from 1990 and 2010 by using the IR-MAD approach[52,57].Only the areas on the Landsat images from 1990 and 2010 that were masked by the changed area were classified with the SVM classifier.The land cover for 2000 was updated with the classification results of the changed area to develop the land cover for 1990 and 2010.

Figure 2 .
Figure 2. Flow diagram of the land cover and land use change (LCLUC) classification.IR-MAD, iteratively re-weighted multivariate alteration detection.

Figure 2 .
Figure 2. Flow diagram of the land cover and land use change (LCLUC) classification.IR-MAD, iteratively re-weighted multivariate alteration detection.

of 21 Figure 3 .
Figure 3. Flow diagram of the land cover change classification.

Figure 3 .
Figure 3. Flow diagram of the land cover change classification.

Figure 4 .
Figure 4. Details of the classifications of the sub-area in the JJT: (a) Landsat 5 TM sub-image at Path/Row 123/032 that was acquired on 7 September 1992; (b) Landsat 7 ETM+ sub-image that was acquired on 2 August 1999; (c) Landsat 5 TM sub-image that was acquired on 22 September 2009; (d-f) classification results for (a-c), respectively; and (g) China's Land-Use/cover Datasets (CLUDs) in 2000.

Figure 4 .
Figure 4. Details of the classifications of the sub-area in the JJT: (a) Landsat 5 TM sub-image at Path/Row 123/032 that was acquired on 7 September 1992; (b) Landsat 7 ETM+ sub-image that was acquired on 2 August 1999; (c) Landsat 5 TM sub-image that was acquired on 22 September 2009; (d-f) classification results for (a-c), respectively; and (g) China's Land-Use/cover Datasets (CLUDs) in 2000.

Figure 5 .
Figure 5. Objects that were segmented by the mean shift method and the refined impervious surface after synthesizing the object and SVM results: (a) ETM+ sub-image at Path 123/Row 032 that was acquired on 21 October 1999; (b) boundaries of objects that were segmented by the mean shift method (yellow lines); (c) impervious surface that was classified by the SVM classifier; (d) refined impervious surface based on the integration of the segmented objects and the SVM classification result; and (e) detailed information regarding the integration.

Figure 5 .
Figure 5. Objects that were segmented by the mean shift method and the refined impervious surface after synthesizing the object and SVM results: (a) ETM+ sub-image at Path 123/Row 032 that was acquired on 21 October 1999; (b) boundaries of objects that were segmented by the mean shift method (yellow lines); (c) impervious surface that was classified by the SVM classifier; (d) refined impervious surface based on the integration of the segmented objects and the SVM classification result; and (e) detailed information regarding the integration.

21 Figure 6 .
Figure 6.Changed area that was determined by comparing two images: (a) ETM+ sub-image at Path 123/Row 032 that was acquired on 2 August 1999; (b) Landsat 5 TM sub-image that was acquired on 22 September 2009; (c) chi-square image from the IR-MAD algorithm; and (d) mask of the changed area, which is shown in yellow against the background image (i.e., Figure 6a).

Figure 6 .
Figure 6.Changed area that was determined by comparing two images: (a) ETM+ sub-image at Path 123/Row 032 that was acquired on 2 August 1999; (b) Landsat 5 TM sub-image that was acquired on 22 September 2009; (c) chi-square image from the IR-MAD algorithm; and (d) mask of the changed area, which is shown in yellow against the background image (i.e., Figure 6a).

Figure 7 .
Figure 7. Classifications of the Beijing-Tianjin-Tangshan (BTT) region for 1990, 2000 and 2010: (a) land cover in 1990; (b) land cover in 2000; and (c) land cover in 2010.The area that is indicated by the black rectangle is shown in detail in Figure 4.

Figure 7 .
Figure 7. Classifications of the Beijing-Tianjin-Tangshan (BTT) region for 1990, 2000 and 2010: (a) land cover in 1990; (b) land cover in 2000; and (c) land cover in 2010.The area that is indicated by the black rectangle is shown in detail in Figure 4.

Figure 8 .
Figure 8. Percentages of different land cover areas for the three periods.

Figure 8 .
Figure 8. Percentages of different land cover areas for the three periods.

Table 1 .
TM/ETM+ images that were used in this study.The acquisition dates comprise the year, month and day.All of the images for 2000 were ETM+ images, whereas those for 1990 and 2010 were TM images.The images that were acquired on 27 August 1993 for 121/33 were used to complement areas that were influenced by cloud cover.

Table 2 .
Adjusted accuracy assessment of the classifications for the three periods."Overall" is the overall accuracy, "Kappa" is the Kappa coefficient, "PA" is the producer's accuracy and "UA" is the user's accuracy.

Table 3 .
Changes in impervious surface after integration based on the sub-image.

Table 3 .
Changes in impervious surface after integration based on the sub-image.

Table 4 .
Error-adjusted areas and proportions of the classes for the three periods.

Table 4 .
Error-adjusted areas and proportions of the classes for the three periods.