Mapping Typical Urban LULC from Landsat Imagery without Training Samples or Self-Defined Parameters

Land use/land cover (LULC) change is one of the most important indicators in understanding the interactions between humans and the environment. Traditionally, when LULC maps are produced yearly, most existing remote-sensing methods have to collect ground reference data annually, as the classifiers have to be trained individually in each corresponding year. This study presented a novel strategy to map LULC classes without training samples or assigning parameters. First of all, several novel indices were carefully selected from the index pool, which were able to highlight certain LULC very well. Following this, a common unsupervised classifier was employed to extract the LULC from the associated index image without assigning thresholds. Finally, a supervised classification was implemented with samples automatically collected from the unsupervised classification outputs. Results illustrated that the proposed method could achieve satisfactory performance, reaching similar accuracies to traditional approaches. Findings of this study demonstrate that the proposed strategy is a simple and effective alternative to mapping urban LULC. With the proposed strategy, the budget and time required for remote-sensing data processing could be reduced dramatically.


Urbanization in Developing Countries
Land use/land cover (LULC) change is one of the most important indicators in understanding the impact of human activities on the environment, which seem unprecedented despite profoundly affecting the Earth's ecological systems [1].Of the LULC changes induced by humans, urbanization manifests itself as the most widespread anthropogenic causes of the loss of arable land, habitat destruction and the decline in natural vegetation cover [2].Global urban areas are rapidly expanding and nowadays, the majority of the world's population lives in cities [2,3].As ecosystems in urban areas are strongly influenced by anthropogenic activities, considerable attention is currently directed towards monitoring changes in urban LULC [4].Such studies are particularly important because the spatial characteristics of LULC are useful for understanding various impacts of human activities on ecological conditions in urban environment.
It has been projected that most of the world's megacities will be in developing countries by 2020, such as China, India, Vietnam and Bangladesh [3].In developing countries, urban expansion is often rapid and unplanned, which can lead to unintended and detrimental consequences [1][2][3][4][5][6][7][8][9][10].Cities are often located on the most productive agricultural lands, so any expansion of built-up areas quickly consumes natural resources, compromising not only food production, but also the provision of ecosystem goods and services that are derived from these landscapes.In developing countries, such as China, India, Bangladesh and Vietnam, state-led industrialization and urban growth policies have often compelled farmers to sell high-yield cropland to developers [2].According to the Ministry of Land and Resources of China, the total cultivated land in Mainland China continually decreased by approximately 5.6 million hectares during last two decades [11].
Furthermore, urbanization also ultimately results in a considerable reduction in other natural lands (e.g., green space, open spaces and water bodies).As the focal point of the urban ecological system, green spaces play a crucial role in cleaning the air, adjusting the microclimate, eliminating noise, beautifying the surroundings and so on.They support the construction of high-quality human settlements, since green spaces act as the "lungs" of the city [3].Furthermore, urbanization often leads to the destruction of local environment, particularly in developing countries.The negative impact of urbanization on the environment include increasing vulnerability to natural hazards, channel-bank and road-surface erosion, habitat destruction, landscape degradation and fragmentation, climate change, species extinction as well as the reduction of net primary productivity [1][2][3][4][5][6][7][8][9][10].
Collectively, the unchecked development of hundreds of cities and towns has raised national and global concerns about energy security, greenhouse gas emissions, environmental changes and major modifications to the natural landscape.In order to mitigate the detrimental effects associated with urban growth on the natural resources, understanding the environment and how to maintain optimal ecosystem functioning, spatial and temporal LULC patterns as well as the factors affecting these changes are considerably important in developing rational economic, social and environmental policies [5].

Remote-Sensing for Analyzing LULC Change
Remote-sensing is a powerful and cost-effective data source for assessing the spatial and temporal dynamics of LULC [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17].Digital archives of remotely sensed data provide an excellent opportunity to study historical land use/cover changes and to relate their spatio-temporal patterns to environmental and human factors.With the rapid development of Earth observation techniques, it has become convenient to obtain a large number of remotely-sensed imagery over a certain area at different times, from hundreds of Earth observation platforms [12].However, this brings challenges to researchers to timely process this large amount of remote-sensing data as well as to rapidly transfer data into information and knowledge [13].Recently, the contradiction between rapid development of data acquisition and the relatively lagging of data processing has become increasingly prominent [14].
In the past decades, researchers are inclined to employ supervised classification methods to obtain LULC maps [18].However, despite the need for substantial expert inputs, these algorithms are not able to account for the variation in imaging dates, image quality, vegetation growth, weather conditions and cloud cover [19].This shortcoming resides in the instantaneous character of remote-sensing images, which only provide snapshots of surface conditions.Classification algorithms and parameters assigned for input data at a certain time are probably not suitable for those at other times, no matter how well they were developed [20].Essentially, samples for training classifiers and parameters have to be collected each time for LULC mapping.National or regional LULC maps are often produced in a yearly basis.In these cases, most existing methods have to collect ground reference data at an annual frequency to train the classifiers [21][22][23].However, collecting reference data and selecting sample sets always require a considerable amount of time and a large budget, which makes it difficult to process a large number of images in a time-sensitive matter.Consequently, monitoring LULC changes across multiple years becomes a hard task, as time series images must be processed.
A few studies have addressed an approach to mapping LULC in multiple years without real-time training samples, namely the classification extension or generalization [24][25][26].This method tries to find sample sets to train classifiers, which is then utilized repeatedly for multiple years without the need for year-to-year reformulation.For example, a phenology-based classification method was used to map corn and soybean in multiple years using training data limited to a single year [19].However, successful examples of classifier extension are rare and mostly occur in natural vegetation environments [19,[25][26][27].
To improve the performance of the classification extension, new sets of input variables with better mapping capability need to be developed.Remote-sensing index has been commonly believed as a simple and effective way to highlight certain LULC types [28,29].For example, vegetation indices (i.e., the Normalized Difference Vegetation Index (NDVI) [30,31] and Enhanced Vegetation Index (EVI) [32]) have been successfully implemented in time-series image analysis [33].As we know, multiple indices for typical urban LULC have been developed, such as the Normalized Difference Water Index (NDWI) [34,35], the Normalized Difference Built-up Index (NDBI) [36] and the Normalized Difference Bare Land Index (NDBaI) [37].In these cases, it is possible to map typical LULC with pre-determined thresholds of these indices.However, it is commonly recognized that assigning a suitable threshold to identify a certain class from an index image is difficult [36,38].A low threshold excludes pixels of this class, while a high value includes too many pixels of other classes.In most studies, repeated experiments and comparisons are in demand to specify the threshold, which requires substantial labor and a considerable amount of time [39].
This study aims to eliminate the need to collect training samples for mapping LULC in multiple years with integrated usage of several indices.The classification strategy is designed as follows.Remote-sensing indices are carefully selected to distinguish specific LULC.With these indices, an unsupervised classifier is employed to build LULC clusters without training samples or assigning thresholds.Following this, a supervised classification is tested with samples randomly collected from the unsupervised clusters.Both are carefully compared with the results from a traditional method.At last, characteristics and potential applications about the proposed strategy are discussed.

Study Site
Wuhan is known as the "Chicago of China" [40,41] due to its critical role in inland transportation.The city's climate is humid subtropical , with plentiful rainfall mainly from May to July [42].As one of the greatest developing cities in China, Wuhan is undergoing the most dramatic expansion in the past decades.Since 2007, the city has had more than 5000 construction fields each year, with many of them lasting more than one year [43].On one hand, the rapid expansion brings great benefits to citizens, such as job opportunities, convenient transportation, beautiful views and modern dwelling environments.On the other hand, it induces serious air pollution (i.e., smog), water pollution and soil loss in the city [44].

Data
Landsat images covering the Wuhan City (Path: 123, Row: 39) were carefully selected for this study (Table 1).The data were downloaded with special consideration for cloud cover, phenology and dryness of ground.In this study, Landsat TIR bands were rescaled to 30 m prior to building an index image.As shown in the subset in Figure 1, many large bare lands are distributed around built-up areas and lakes, illustrating on-going urban expansion.The basic maps, climate and environment data of the city were collected from Wuhan Municipal Environmental Protection Bureau.Related social, commercial and transportation data came from Wuhan Statistical Yearbook.
An independent set of very high resolution (VHR) imagery (IKONOS) covering the study site in 2000, 2004, 2007, 2009, 2011 and 2013 were collected and served as validation sources.The manually interpreted bare lands from the VHR images were assumed to represent ground truth.The stratified randomized method was adopted to generate two hundred samples for both the bare land and "the other".
Atmospheric correction of all Landsat images in Table 1 was implemented by the FLAASH model in ENVI 5.1, using the Landsat sensor response functions, metadata of the Landsat images and meteorological data at the Wuhan station.Digital numbers at each pixel were converted to surface reflectance through this process.

Indices Related to Typical LULC
In this study, urban LULC lands are classified into five typical types, namely built-up, bare land, water body, agricultural land and forest.In this present study, agricultural land is actually a mixture of crop fields, grasslands and other herbaceous covers.For each LULC type, a large number of remote-sensing indices have been published [30,32].This study summarized these indices in Table 2.The basic maps, climate and environment data of the city were collected from Wuhan Municipal Environmental Protection Bureau.Related social, commercial and transportation data came from Wuhan Statistical Yearbook.
An independent set of very high resolution (VHR) imagery (IKONOS) covering the study site in 2000, 2004, 2007, 2009, 2011 and 2013 were collected and served as validation sources.The manually interpreted bare lands from the VHR images were assumed to represent ground truth.The stratified randomized method was adopted to generate two hundred samples for both the bare land and "the other".
Atmospheric correction of all Landsat images in Table 1 was implemented by the FLAASH model in ENVI 5.1, using the Landsat sensor response functions, metadata of the Landsat images and meteorological data at the Wuhan station.Digital numbers at each pixel were converted to surface reflectance through this process.

Indices Related to Typical LULC
In this study, urban LULC lands are classified into five typical types, namely built-up, bare land, water body, agricultural land and forest.In this present study, agricultural land is actually a mixture of crop fields, grasslands and other herbaceous covers.For each LULC type, a large number of remote-sensing indices have been published [30,32].This study summarized these indices in Table 2. Vegetation indices are usually the ratios or linear combinations of reflectance measurements in different spectral bands, especially the visible and near-infrared bands, which are used to enhance spectral differences between green vegetation and other LULCs.A considerable number of vegetation indices have been proposed, ranging from very simple to very complex band combinations, including the Band Ratio (RATIO), Band Difference (DVI), Band Sum (SUM), the Vegetation Index (VI), and the Transformed Vegetation Index (TVI), Soil-Adjusted Vegetation Index (SAVI) [49], enhanced vegetation index (EVI), Normalized Difference Vegetation Index (NDVI) and so on.
Water indices are developed to delineate open water features.For example, the Normalized Difference Water Index (NDWI) maximizes reflectance of water at green light channel and minimizes its low reflectance at NIR [35].It also takes advantage of high reflectance of vegetation and soil features at NIR.A previous study [34] proposed another NDWI to study water content in vegetation, using a near infrared channel (centered approximately at 0.86 µm) and a mid-infrared channel (at 1.24 µm).This index can effectively reflect the water stress of canopy.Another study [28] modified the NDWI by using a middle infrared (MIR) band, such as TM5, to substitute the NIR band in it.In the modified NDWI (MNDWI), noises, such as built-up pixels, are successfully eliminated from the enhanced water features.
Built-up lands are dominated by impervious surfaces.Conversion from nature lands to impervious lands often results in serious negative outcomes, such as heat islands, waterlogging and photo-chemical smog.Many researchers have made use of remote-sensing imagery to extract built-up lands.One study [36] proposed the Normalized Difference Built-up Index (NDBI) using TM4 and TM5 and applied it in extracting urban areas of Nanjing City of China.Another study [46] proposed the Urban Index (UI) using Landsat TM7 and TM4, exploiting an observed inverse relationship between the brightness of urban areas in near-infrared (0.76-0.90 µm) and mid-infrared (2.08-2.35µm).The Index-based Built-Up Index (IBI) [47] made use of three thematic indices rather than original spectral bands to construct a new index.The subtraction of the SAVI band and the MNDWI band from the NDBI band results in positive values for built-up pixels.
Due to similar spectral responses of bare lands and built-ups in addition to the presence of mixed pixels in urban area, most built-up indices fail to delineate the two LULC types in a clear manner [38], although some studies have used grey level co-occurrence matrices to extract urban areas with good results [50].In a developing city, bare lands should not be misclassified, owing to its large coverage, unique physical characters and environmental impacts.Recently, one study [37] proposed the Normalized Difference Bare Land Index (NDBaI) to distinguish the bare land, with Landsat band 5 (SWIR1) and band 6 (TIR).Another study [35] presented the Enhanced Built-up and Bare Land index (EBBI) to map the built-up and bare land at the same time, using the NIR, SWIR1, and TIR bands (Landsat ETM+ bands 4, 5 and 6).One study [48] proposed a Normalized Bare Land Index (NBLI) to distinguish bare land areas, using the Red and Blue wavelength, which correspond to band 3 and band 6 in Landsat TM.
As listed in Table 2, all these indices to extract vegetation, water, built-up and bare lands were compared in this study.

Comparative Analysis
In order to find the optimal index to distinguish a specific LULC type, the capacity of all indices in Table 2 were comparatively analyzed using hundreds of samples selected from index images (Figure 2).
Remote Sens. 2017, 9, 700 6 of 23 (EBBI) to map the built-up and bare land at the same time, using the NIR, SWIR1, and TIR bands (Landsat ETM+ bands 4, 5 and 6).One study [48] proposed a Normalized Bare Land Index (NBLI) to distinguish bare land areas, using the Red and Blue wavelength, which correspond to band 3 and band 6 in Landsat TM.
As listed in Table 2, all these indices to extract vegetation, water, built-up and bare lands were compared in this study.

Comparative Analysis
In order to find the optimal index to distinguish a specific LULC type, the capacity of all indices in Table 2 were comparatively analyzed using hundreds of samples selected from index images (Figure 2).In the Water index category in Figure 2, all three indices are able to enhance water features, while suppressing vegetation and soil types.However, built-up lands were found to have relatively high values in the NDWI1 image.Studies have shown that the extracted water was often mixed up with built-up land noises [41].In the NDWI2 image, water is easily confused with vegetation.They have similar spectral differences between the NIR band (TM4) and MIR band (TM5).Thus, the calculation of NDWI2 produces high values for both vegetation and water.In MNDWI equation, the built-up and vegetation pixels have lower values while water has higher values.Therefore, MNDWI maximally delineated water from other classes.
In the Vegetation index category, vegetation features are enhanced, while others are suppressed in both NDVI and DVI image.However, agricultural lands and forest have similar values in these two indices because of their spectral similarity.Thus, traditional vegetation indices cannot be used to distinguish the forest from agricultural land.In the Water index category in Figure 2, all three indices are able to enhance water features, while suppressing vegetation and soil types.However, built-up lands were found to have relatively high values in the NDWI1 image.Studies have shown that the extracted water was often mixed up with built-up land noises [41].In the NDWI2 image, water is easily confused with vegetation.They have similar spectral differences between the NIR band (TM4) and MIR band (TM5).Thus, the calculation of NDWI2 produces high values for both vegetation and water.In MNDWI equation, the built-up and vegetation pixels have lower values while water has higher values.Therefore, MNDWI maximally delineated water from other classes.
In the Vegetation index category, vegetation features are enhanced, while others are suppressed in both NDVI and DVI image.However, agricultural lands and forest have similar values in these two indices because of their spectral similarity.Thus, traditional vegetation indices cannot be used to distinguish the forest from agricultural land.
In the Built-up index category, the built-up pixels are mixed up with bare land noises as both have similar values from the Green (TM2) to MIR band (TM5).It is impossible to clearly distinguish built-up from bare lands using those indices.Furthermore, vegetation (forest and agriculture) has values closer to the built-up land in NDBI and IBI, which may be related to higher MIR reflectance from drier vegetation.In the UI image, the distance between the built-up land and vegetation becomes further.
In the bare land index category, all indices are able to enhance the bare land features.In NDBaI, bare lands are probably mixed up with vegetation and built-up noises.In EBBI, the error range of built-up noises also overlaps with bare land, although the distance between them becomes slightly larger than in NDBaI.In NBLI, bare land is separated clearly from all others.In particular, the distance between bare land and built-up areas becomes far enough to be easily separated.This indicates a promising tool to distinguish bare land from built-up areas in a clear manner.The unique problem is that the error range of water may partly cover that of bare land.To avoid this problem, water needs to be removed ahead of time.
It is also noticed that the difference between agricultural land and forest is clear in NBLI, as the index is able to highlight soil composition at a pixel level.As the tree canopies are much denser than crops, tree pixels have less soil effects than crop pixels.Therefore, it is possible to separate forest from agricultural land with the NBLI.To highlight the forest, which has the lowest value in NBLI, an inverse NBLI index is generated as (TIR − Red)/(TIR + Red).

Unsupervised Classification with Selected Indices
The comparison analysis above indicates that some indices have the ability to enhance a certain LULC, especially MNDWI, NBLI and UI.Due to the apparent differences of LULC types in these indices, it is possible to classify these types without training samples or thresholds.An unsupervised classifier, the K-means approach, is employed to cluster the index images with default settings.In each index image, the class type with the highest values is automatically assigned as the target.As shown in Figure 3, a stepwise mapping method is proposed, which includes the following steps.
Step 1: Constructing the MNDWI, NBLI and UI image; Step 2: Extracting water from the MNWI image with the K-means classifier; Step 3: Removing water in the NBLI image, then extracting bare land with the K-means classifier; Step 4: Removing water and bare land in the UI image, then extracting built-up with the K-means classifier; Step 5: Removing above classes in the inverse NBLI image, then extracting forest with K-means classifier; Step 6: Labeling left pixels as agricultural lands; Step 7: Composite all layers into an initial LULC class map, and building the confusion matrix for accuracy assessment.
Specifically, when a land cover was removed from an index image, a negative number was given to the corresponding areas.Following this, the unsupervised classifier was implemented on this image.The class with the highest value was extracted, mapped and given fixed class number, such as 1 for built-up areas, 2 for forest and so on.After this, this thematic map was employed to remove land cover from another index image for extracting another class.At last, outputs in the previous steps were stacked and layered into the final LULC map.

Supervised Classification with Samples Selected Automatically
In the initial LULC map above, pixels with high values in an index image have a high possibility of representing the related LULC.Under this hypothesis, an automatic method is presented to select representative samples from the initial map from the unsupervised classification.
Step 1: Select all pixels of a certain LULC from the initial map and calculate its basic statistics of the related index (i.e., the mean and standard deviation ).
Step 2: With a stratified randomized sampling method, collect sample pixels for each LULC with its related index values located in the range from [ + ] to [ + 3 ].The range is assigned in this study to filter out extreme points (noises) but remain sufficient training samples.If a researcher is inclined to select "pure samples" or "end members", a small range from [ + 2 ] to [ + 3 ] is suggested.
Step 3: Implement supervised classification with a typical classifier to map urban LULC.
Step 4: Build the confusion matrix for accuracy assessment with validation samples.The traditional confusion (error) matrix for a single-date classification was employed in this study.From the confusion matrix, the overall accuracy was calculated by dividing the sum of the entries along the diagonal by the total number of validation samples.The Kappa coefficient of agreement was also derived from the confusion matrix.

Supervised Classification with Samples Selected Automatically
In the initial LULC map above, pixels with high values in an index image have a high possibility of representing the related LULC.Under this hypothesis, an automatic method is presented to select representative samples from the initial map from the unsupervised classification.
Step 1: Select all pixels of a certain LULC from the initial map and calculate its basic statistics of the related index (i.e., the mean E and standard deviation σ).
Step 2: With a stratified randomized sampling method, collect sample pixels for each LULC with its related index values located in the range from [E + σ] to [E + 3σ].The range is assigned in this study to filter out extreme points (noises) but remain sufficient training samples.If a researcher is inclined to select "pure samples" or "end members", a small range from [E + 2σ] to [E + 3σ] is suggested.
Step 3: Implement supervised classification with a typical classifier to map urban LULC.
Step 4: Build the confusion matrix for accuracy assessment with validation samples.The traditional confusion (error) matrix for a single-date classification was employed in this study.From the confusion matrix, the overall accuracy was calculated by dividing the sum of the entries along the diagonal by the total number of validation samples.The Kappa coefficient of agreement was also derived from the confusion matrix.
The results in Figure 4 should be compared with the true color image in Figure 1 to analyze the performance of indices on enhancing certain LULC.In the MNDWI image (Figure 4a), water is evidently enhanced and all other types are compressed effectively.In the NDWI1 image (Figure 4b), gray tones in built-up and bare lands are close to water.In NDWI2 image (Figure 4c), the values of vegetation lands (agriculture and forest) also get close to water.Therefore, it becomes difficulty to extract water accurately from the latter two images.
From Figure 4d,e, agricultural land and forest have similar values in NDVI and DVI images.It can be concluded that typical vegetation index cannot be used to separate agricultural lands from forest.In the inverse NBLI image (Figure 4k), the differences between agricultural land and forest are much clearer.
In Figure 4f-h, bare land and built-up areas share the same gray range.In this case, it is impossible to separate each other, no matter which approaches we take.That is to say, these indices cannot be directly employed to extract built-up areas in developing cities, as the bare land cannot be ignored there.It is also noticed that both built-up and bare land get further away from other classes in the UI image (Figure 4h).Therefore, it is possible to extract built-up areas if bare land could be masked out of the UI image.
In the NDBaI image (Figure 4i), the tonal difference between bare land and built-up is clear.However, some other classes, such as agricultural land, get very close to bare land, which may induce serious misclassification.In the NBLI image (Figure 4k), the difference between bare land and other classes is very clear.However, some water bodies have similar values with bare land, like the Yangzi River, due to the suspended matter in water.
Remote Sens. 2017, 9, 700 9 of 23 The results in Figure 4 should be compared with the true color image in Figure 1 to analyze the performance of indices on enhancing certain LULC.In the MNDWI image (Figure 4a), water is evidently enhanced and all other types are compressed effectively.In the NDWI1 image (Figure 4b), gray tones in built-up and bare lands are close to water.In NDWI2 image (Figure 4c), the values of vegetation lands (agriculture and forest) also get close to water.Therefore, it becomes difficulty to extract water accurately from the latter two images.
From Figure 4d,e, agricultural land and forest have similar values in NDVI and DVI images.It can be concluded that typical vegetation index cannot be used to separate agricultural lands from forest.In the inverse NBLI image (Figure 4k), the differences between agricultural land and forest are much clearer.
In Figure 4f-h, bare land and built-up areas share the same gray range.In this case, it is impossible to separate each other, no matter which approaches we take.That is to say, these indices cannot be directly employed to extract built-up areas in developing cities, as the bare land cannot be ignored there.It is also noticed that both built-up and bare land get further away from other classes in the UI image (Figure 4h).Therefore, it is possible to extract built-up areas if bare land could be masked out of the UI image.
In the NDBaI image (Figure 4i), the tonal difference between bare land and built-up is clear.However, some other classes, such as agricultural land, get very close to bare land, which may induce serious misclassification.In the NBLI image (Figure 4k), the difference between bare land and other classes is very clear.However, some water bodies have similar values with bare land, like the Yangzi River, due to the suspended matter in water.

Unsupervised Classification
In this experiment, three selected indices (i.e., MNDWI, NBLI, UI and inversed NBLI) and the K-means classifier with default settings (i.e., 4 classes with 100 iterations at maximum) were adopted to extract initial maps of typical urban LULC automatically.The initial results are shown in Figure 5 and Table 2.Only one targeted class was extracted in each index image.All non-target classes were combined and labeled as "Other".

Unsupervised Classification
In this experiment, three selected indices (i.e., MNDWI, NBLI, UI and inversed NBLI) and the K-means classifier with default settings (i.e., 4 classes with 100 iterations at maximum) were adopted to extract initial maps of typical urban LULC automatically.The initial results are shown in Figure 5 and Table 2.Only one targeted class was extracted in each index image.All non-target classes were combined and labeled as "Other".
Figure 5 illustrates that it is feasible to extract typical urban LULC types with the proposed workflow as stated in Section 3.3.Aside from water, bare land and built-up areas, the method could also separate forest from agricultural areas effectively without assigning samples or parameters.Comparing Figure 5f to Figure 5b, many built-up areas are misclassified as bare land due to the Ostu result.This proves that binary classifiers, such as the Otsu, are unable to find the right threshold for a certain class, as there were often more than two classes in the index image.Figure 5 illustrates that it is feasible to extract typical urban LULC types with the proposed workflow as stated in Section 3.3.Aside from water, bare land and built-up areas, the method could also separate forest from agricultural areas effectively without assigning samples or parameters.Comparing Figure 5f to Figure 5b, many built-up areas are misclassified as bare land due to the Ostu result.This proves that binary classifiers, such as the Otsu, are unable to find the right threshold for a certain class, as there were often more than two classes in the index image.
In Table 3, the overall accuracy and kappa coefficient achieve 88.95% and 0.8619 respectively.The overall indicators similar to those produced by traditional methods, and satisfy the requirements of most applications.For most classes, their producer's accuracy and user's accuracy are higher than 90%.Due to spectral similarity between agriculture and forest, a lower user accuracy of agriculture (69.65%) and a lower producer's accuracy of forest (66.02%) were resulted.These accuracies are acceptable in this study since classification in vegetative lands is not the focus.In Table 3, the overall accuracy and kappa coefficient achieve 88.95% and 0.8619 respectively.The overall indicators similar to those produced by traditional methods, and satisfy the requirements of most applications.For most classes, their producer's accuracy and user's accuracy are higher than 90%.Due to spectral similarity between agriculture and forest, a lower user accuracy of agriculture (69.65%) and a lower producer's accuracy of forest (66.02%) were resulted.These accuracies are acceptable in this study since classification in vegetative lands is not the focus.

Supervised Classification
Relying on the initial map of unsupervised classification in Figure 5e, training samples were statistically identified and the support vector machine (SVM) classifier was adopted for advanced classification of urban LULC types.The results are shown in Figure 6.

Supervised Classification
Relying on the initial map of unsupervised classification in Figure 5e, training samples were statistically identified and the support vector machine (SVM) classifier was adopted for advanced classification of urban LULC types.The results are shown in Figure 6.It is evident that the ASC and TSC results match those in Figure 4e very well, which indicates that the three methods have similar overall performance.Particularly, both supervised methods better distinguished the forest from agricultural land.For example, the Luojia hill (a famous scenic site of Wuhan University, marked by the black ellipse in the figure) becomes much clearer and more complete in Figure 6b.
As shown in Table 4, the ASC results also reach satisfied overall accuracy and kappa coefficient, similar to that of the unsupervised classification).Specifically, the misclassification between built-up and bare land becomes serious in ASC result, as the producer's accuracy of bare land and the user's It is evident that the ASC and TSC results match those in Figure 4e very well, which indicates that the three methods have similar overall performance.Particularly, both supervised methods better distinguished the forest from agricultural land.For example, the Luojia hill (a famous scenic site of Wuhan University, marked by the black ellipse in the figure) becomes much clearer and more complete in Figure 6b.
As shown in Table 4, the ASC results also reach satisfied overall accuracy and kappa coefficient, similar to that of the unsupervised classification).Specifically, the misclassification between built-up and bare land becomes serious in ASC result, as the producer's accuracy of bare land and the user's accuracy of built-up become lower (75.25% and 78.17% respectively).On the other hand, the difference between forest and agriculture becomes much clearer in the result, as all related indicators are higher than 80%.Table 5 shows that the TSC has very good performance.The overall accuracy and kappa efficiency achieve 94.09% and 0.9261 respectively, and all detailed indicators are close to or higher than 90%.It is noted that the accuracy of both UC and ASC productions are very close to that of TSC result.That is to say, the proposed automatic approach could produce results with satisfied accuracy, so that the method could be considered as an alternative or supplemental option to existing methods.

Time Series Image Analysis
With the proposed method, the LULC maps of Wuhan city in 2000,2004,2007,2009,2011 and 2013 are generated to examine the LULC changes since 2000 (Figure 7).The post-classification comparison of change detection was carried out using ENVI 5.1, producing a summary of the major LULC conversions, namely 'from-to' information, which occurred during the study period.
Spatial patterns of LULC changes in the Wuhan area for 2000, 2004, 2007, 2011 and 2013 are shown in Figure 7.The city has expanded dramatically since 2000 in both built-up areas and bare lands.Analysis of the LULC changes in Table 6 revealed an increase in these land types over the study period.Built-up areas increased by 71.52 km 2 between 2000 and 2004; 120.12 km 2   2 between 2011 and 2013.The net increase in bare land areas over the study period was 167.91 km 2 , which is an average of more than 12.9 km 2 /year.When compared with other cities in developing counties, such as Dhaka City in Bangladesh, the rate of the built-up expansion in Dhaka City was 4.47 km 2 /year, while that of bare land was 4.24 km 2 /year between 2000 and 2011 [10].The nature of urban expansion is generally related to demographic change, economic growth, topography, land use and transportation [1].Close examination of the change detection statistics revealed that approximately 648.80 km 2 of the built-up area and bare land in Wuhan were previously either agricultural areas between 2000 and 2013, 34.39 km 2 of them were from water bodies and 37.51 km 2 of them were from forest.
Remote Sens. 2017, 9, 700 15 of 23 previously either agricultural areas between 2000 and 2013, 34.39 km 2 of them were from water bodies and 37.51 km 2 of them were from forest.Generally, two national economic zones were observed to have promoted built-up and bare land growth.In the southwest, the main development is the national Wuhan Economic Development Zone dominated by the vehicle production industry.In the east, the main development area is the national East Lake High Tech Development Zone known as the "China's optic valley".The two zones are engines driving the city's economic development.Large areas of agricultural lands have been exploited for constructions, such as that in Hanoi [8].
The analysis also revealed that the area occupied by water bodies decreased by 2.75%, forest by 8.18% and agricultural land by 12.52% between 2000 and 2013.Significant increases in built-up and bare land were observed in this period.Specifically, built-up areas increased in size by 204.45% and bare land by 40.96% in the period from 2000 to 2013.The decline of vegetation and water bodies was clearly due to intensification of urban development in the Wuhan, particularly due to the establishment of national development zones.

Another Example
Guangzhou is currently the third highest populated city in mainland China, behind Beijing and Shanghai.It is located within the Pearl River Delta at the downstream of Pearl River, about 120 km (75 mi) north-northwest of Hong Kong and 145 km (90 mi) north of Macau.In 2015 the city's administrative area was estimated to have a population of 13,501,100 and belongs to one of the most populous metropolitan agglomerations on Earth.In 2008, the migrant population from other provinces of China in Guangzhou was 40 percent of the city's total population.For the three consecutive years 2013-2015, Forbes ranked Guangzhou as the best commercial city on the Chinese  Generally, two national economic zones were observed to have promoted built-up and bare land growth.In the southwest, the main development is the national Wuhan Economic Development Zone dominated by the vehicle production industry.In the east, the main development area is the national East Lake High Tech Development Zone known as the "China's optic valley".The two zones are engines driving the city's economic development.Large areas of agricultural lands have been exploited for constructions, such as that in Hanoi [8].
The analysis also revealed that the area occupied by water bodies decreased by 2.75%, forest by 8.18% and agricultural land by 12.52% between 2000 and 2013.Significant increases in built-up and bare land were observed in this period.Specifically, built-up areas increased in size by 204.45% and bare land by 40.96% in the period from 2000 to 2013.The decline of vegetation and water bodies was clearly due to intensification of urban development in the Wuhan, particularly due to the establishment of national development zones.

Another Example
Guangzhou is currently the third highest populated city in mainland China, behind Beijing and Shanghai.It is located within the Pearl River Delta at the downstream of Pearl River, about 120 km (75 mi) north-northwest of Hong Kong and 145 km (90 mi) north of Macau.In 2015 the city's administrative area was estimated to have a population of 13,501,100 and belongs to one of the most populous metropolitan agglomerations on Earth.In 2008, the migrant population from other provinces of China in Guangzhou was 40 percent of the city's total population.For the three consecutive years 2013-2015, Forbes ranked Guangzhou as the best commercial city on the Chinese mainland.Guangzhou has a warm, monsoon-influenced, humid subtropical climate.The primary LULCs in the city include built-up, forest, water and bare land.Our test is implemented with a Landsat image acquired on 2 November 2009.The results are shown in Figure 8.   Figure 8 reveals that the proposed method could be successfully implemented in a different region.The differences between water, bare land, built-up areas, forest and agricultural land are clear in the index images.The classification results match the visual interpretation very well, and the overall accuracy and kappa coefficient achieve 89.2% and 0.8734 respectively.Only few great bare lands were observed within the city core, illustrating that Guangzhou has been well developed in 2009.Contrarily, Wuhan is still under rapid constructing.

Characteristics of the Automatic Method
Remote sensing index has been proved as a simple, reliable and effective way to highlight and extract certain LULC.It is a much faster and cheaper approach to map LULC, compared with various classification methods.However, for some special classes, it may be hard to find a suitable index.For example, it is hard to distinguish the forest from crop land with vegetation indices.In this study, several indices were carefully selected from hundreds of published to map typical urban LULC, including MNDWI, NBLI and UI, for highlighting the water, bare land and built-up respectively.Besides, the inverse NBLI was adopted to separate forest from agricultural land.It was thought to be rational and effective to some extent, as the index could reflect the difference in soil percentage of a pixel.Since the tree canopy are much denser than crops', pixels belonging to trees should have lower soil percentage than those belonging to crops.Then, an attempt towards classifying typical urban LULC automatically without assigning training samples or thresholds with proposed indices, were implemented.Experiment results have verified the feasibility and validity of the proposed method.This study points out a promising way to deal with remote sensing big data, especially on mapping urban LULC with Landsat time series images.
As binary classifiers, like the Otsu, are unable to find the right threshold for a certain class, the common unsupervised classifier, K-means, was adopted to define the threshold for extracting specific class from the index image, with default settings.Experiment results prove the method is very effective to extract a certain LULC from related index image.This proves the proposed method is able to obtain satisfied results without special settings or assumptions, as the K-means is a simple and common classifier.It is possible that better results could be obtained if other unsupervised classifiers are adopted in the process.Tests on other unsupervised classifiers may achieve better results.
In the unsupervised classification result, samples for training supervised classifier were randomly selected from the set of pixels having higher value in the index image.It was believed that those pixels could represent the typical samples of related LULC.In order to assign the index value range for defining the set of candidates, tests on different set of pixels were carried out and analyzed.It was found that using samples with very high index value in the supervised classification often leaded to lower accuracy, even if extreme pixels had been filtered.So those samples should not be thought as pure pixels.At last, a suitable range was recommended and stated in Section 3.4.Results illustrated that the proposed methods could achieve satisfied accuracy, similar to the traditional approach.The supervised method could distinguish the forest from agricultural land more clearly.However, its performance became worse when dealing with other classes.Generally, the supervised classification did not improve the accuracy dramatically, so the unsupervised method is recommended for practical implementation.
A problem will appear when dealing with detailed subclasses other than the five classes, that how to find suitable indices for them.Generally, remote sensing indices were mainly developed for typical and evident classes in median or low resolution images, while few indices were designed for subclasses, like grass, corn, soybean, wheat, etc.It is possible to detect insignificant differences between subclasses, by combining multiple indices.For example, water bodies extracted from the MNDWI image, can be further divided according to the percentage of suspended solids in them, which could be displayed in NBLI image to some extent.More analysis and testes on the spectral characters of subclasses, are necessary for this.
Although the method has been successfully applied in the city of Wuhan and Guangzhou, and dealing with time series images, its adaptability should also be further analyzed and tested.Conventionally, classification methods often depends on the local landscape being classified, and landscape really change a lot depending on the part of the globe.It is also noticed that, some special situations may influence the classification accuracy of the proposed method.For example, a follow agricultural land is probably classified as a bare land, and a flooded land may be classified as a lake.To avoid this error, special consideration for cloud, phenology and dryness of ground should be paid.In addition, more tests are needed to verify and improve the adaptability of the proposed method.

Analyses on Urbanization of Wuhan
With the proposed method and time series images, LULC changes of Wuhan have been mapped.Urbanization process of the city was monitored and analyzed from the maps, since urban growth is crucial for local environmental and ecological management.The city of Wuhan has undergone one of the most rapid economic, industrial and urban development in its history, similar to other cities in China, India, Bangladesh, Vietnam and other developing counties [1][2][3][4][5][6][7][8][9][10][51][52][53][54][55][56][57].The study revealed that the urban expansion in Wuhan has been rapid, with the majority of built-up and bare land acquired by converting areas that were previously agricultural land, forest and water bodies.This suggests the existence of increased pressure on natural resources in the city to meet the increasing demand for urban land.
A related study disclosed that LULC changes and urban expansion of the city are governed by a combination of geographical, environmental and socio-economic factors.Rapid growth of the economy is the primary cause for rapid urbanization, as economic opportunities in the city are fueling rural to urban migration on a massive scale [51].Census data indicated that the urban GDP of Wuhan was only 17.76 billion dollars in 2000, which increased to 133.17 billion dollars in 2013, with an annual growth rate of 49.99% [42,43].During the period, the urban population of Wuhan was 7.58 million in 2000, which slightly increased to 8.95 million in 2013, with an annual growth rate of 1.39% [42,43].In Dhaka, Bangladesh, population growth is the main cause for rapid urbanization, as its annual average population growth was 5% during the same period [1].
The rapid and poorly managed expansion, together with the surging popularity of private vehicles, has led to a plague of "city diseases", such as flooding and water logging, air pollution, water pollution, noise, solid waste in addition to ecosystem and environment changes [44].This process is similar that which has occurred in Ho Chi Minh [2], Dhaka [1,[3][4][5][6]10], Delhi [7], Hanoi [8] and Ansan [9].Numerous bare land areas emerging along with urbanization processes have cast irreversible impacts on the urban environment, such as air pollution and soil loss [51][52][53][54][55][56][57].In the plum rain season, flash flooding hazards have become a serious issue and have caused human death and damages to urban infrastructure in the city [52], similar to that in Dhaka [3].One study [53] investigated the consequence of LULC change on thermal environment.This study found that built-up and bare land types consistently exhibited the highest mean LST (Land Surface Temperature), while water-bodies consistently exhibited the lowest temperatures each year and did not significantly change throughout the study period.These findings reflect the results of similar studies in other developing cities.For instance, one study [6] found the increase in unmanaged urbanization in Dhaka and its immediate surroundings has led to a continuous increase in LST as vegetation and floodplains are converted into either bare land or built-up surfaces.
As indicated in this study, urbanization requires more and more built-up areas for housing, business and transport networks, which is generally being met through the development of natural lands (e.g., agriculture lands, forests, water bodies and so on).This ultimately results in a considerable reduction in the open and green areas of that region.Related studies have reported that the increasing level of urbanization was leading to landscape fragmentation in Wuhan [54,55].This resembles the findings in previous studies [4,5] that the rapid human activity in Dhaka, Bangladesh has leaded to numerous unconnected small patches of vegetation and cultivated land, greater isolation as well as higher percentage of edge areas in patches.
Rapid urbanization, along with improper land use, manufacturing industries and large number of vehicles has probably resulted in serious air pollution, water pollution and solid waste [56].According to the report from the local government [57], 160 days in 2013 had good air quality, which only accounts for 43.8% of the year.The annual average concentration of PM2.5 was 94 µg/m 3 , which is 1.7 times higher than the criterion.Only 51.5% of the year had a PM2.5 lower than the criterion.Investigation of water quality disclosed that 23.3% monitored stations of 11 major streams, 14.6% of 89 major lakes and 41.1% of 73 other lakes could not meet the lowest water quality criterion (the Chinese V type) [57].Related monitoring showed the average noise on road at night in downtown was 60 decibels, which is labelled as officially bad [57].Furthermore, the general industrial solid waste was 13,815,500 tons in 2013 [57].Similarly, the report about another developing city, Dhaka Metropolitan Area of Bangladesh, revealed that the increasing population, triggered by rural-urban migration, imposed enormous pressure on the meagre resource bases and resulted in serious environmental degradation of the city [3].

Conclusions
Most developing cities in China and other countries have undergone one of the most rapid urban development in human history.The rapid urbanization in developing cities has also led to an increase in "city diseases", natural hazards and environmental changes.Traditionally, researchers are inclined to map LULC with supervised classification methods, which always needs a considerable amount of time and a large budget.Using remote-sensing indices has been proved as a simple, reliable and effective way.However, assigning a suitable threshold is always difficult.This study presented a novel strategy to map urban typical LULC without training samples or specifying parameters.First of all, several novel indices were carefully selected from the index pool.Following this, an unsupervised classifier was employed to extract LULC from the index image without assigning thresholds.After this, an automatic supervised classification was also tested with samples randomly collected from previous result.
The overall accuracy and kappa coefficient of the unsupervised classification results achieve 88.95% and 0.8619 respectively, which are not much lower than those of traditional methods.Specifically, most producer's accuracies and user's accuracies are higher than 90%, except the user accuracy of agriculture (69.65%) and the producer's accuracy of forest (66.02%).The supervised classification also obtains satisfied overall accuracy and kappa coefficient.Here, all indicators related to forest and agriculture become higher than 80%, while the producer's accuracy of bare land and the user's accuracy of built-up become lower (75.25% and 78.17% respectively).
The proposed strategy can save considerable time and labor forces for mapping LULC from remote-sensing images in multiple years.This will help researchers obtain more findings on the spatial-temporal patterns of LULC changes as well as the interactions between human and surroundings.Further studies will be carried out to find new indices for extracting some subclasses.The study also suggests that it is possible to detect insignificant differences between subclasses by combining multiple indices.This study points out that it is possible to map LULC automatically with several suitable indices and a cluster method, although further analysis on indices and clusters are needed.

Figure 1 .
Figure 1.Study site: the upper left panel shows the location of Wuhan in China; the right panel illustrates the city boundary and an example Landsat image acquired on 31 July 2007; and the bottom left panel shows the detail of an image subset.

Figure 1 .
Figure 1.Study site: the upper left panel shows the location of Wuhan in China; the right panel illustrates the city boundary and an example Landsat image acquired on 31 July 2007; and the bottom left panel shows the detail of an image subset.

Figure 2 .
Figure 2. Statistics of samples of typical LULC types.All indices are grouped into four categories: Water index, Vegetation index, built-up index and Bare land index.All digital numbers were stretched into [0, 255].Marks represent the average value of each type, while the bar in each direction indicates one standard deviation from the average.

Figure 2 .
Figure 2. Statistics of samples of typical LULC types.All indices are grouped into four categories: Water index, Vegetation index, built-up index and Bare land index.All digital numbers were stretched into [0, 255].Marks represent the average value of each type, while the bar in each direction indicates one standard deviation from the average.

Figure 3 .
Figure 3.The flowchart of the unsupervised classification method.

Figure 3 .
Figure 3.The flowchart of the unsupervised classification method.

Figure 5 .
Figure 5. Unsupervised classification results: (a) MNDWI and extracted water in the subset; (b) watermasked NBLI and extracted bare land in the subset; (c) masked UI and extracted built-up in the subset; (d) masked inverse NBLI and the extracted forest and agriculture in the subset; (e) the combination of all LULC; and (f) comparison test, partial masked NBLI and extracted bare land with Ostu.

Figure 5 .
Figure 5. Unsupervised classification results: (a) MNDWI and extracted water in the subset; (b) water-masked NBLI and extracted bare land in the subset; (c) masked UI and extracted built-up in the subset; (d) masked inverse NBLI and the extracted forest and agriculture in the subset; (e) the combination of all LULC; and (f) comparison test, partial masked NBLI and extracted bare land with Ostu.
Remote Sens. 2017, 9, 700 17 of 23 mainland.Guangzhou has a warm, monsoon-influenced, humid subtropical climate.The primary LULCs in the city include built-up, forest, water and bare land.Our test is implemented with a Landsat image acquired on 2 November 2009.The results are shown in Figure 8.

Figure 8 .
Figure 8.An example test in Guangzhou City.(a) The Landsat image of the city on 2 November 2009; (b) the MNDWI image in psudo-color; (c) the NBLI image in psudo-color; (d) the UI image in psudocolor; (e) the inverse NBLI image in psudo-color; and (f) the unsupervised classification result.

Figure 8 .
Figure 8.An example test in Guangzhou City.(a) The Landsat image of the city on 2 November 2009; (b) the MNDWI image in psudo-color; (c) the NBLI image in psudo-color; (d) the UI image in psudo-color; (e) the inverse NBLI image in psudo-color; and (f) the unsupervised classification result.

Table 1 .
Landsat images used in this study.

Table 3 .
Error matrix of the unsupervised classification.

Table 3 .
Error matrix of the unsupervised classification.

Table 4 .
Error matrix of the automatic supervised classification (ASC).

Table 5 .
Error matrix of the traditional supervised classification (TSC).
between 2004 and 2007 km 2 between 2007 and 2009; 128.03 km 2 between 2009 and 2011; and 110.35 km 2 between 2011 and 2013.The net increase in built-up areas over the study period was 552.69 km 2 , which is an average of more than 42 km 2 /year.Similarly, bare land areas increased by 78.36 km 2 between 2000 and 2004; 35.89 km 2 between 2004 and 2007; 16.33 km 2 between 2007 and 2009; 15.41 km 2 between 2009 and 2011; and 21.92 km