Spatiotemporal analysis of urban environment based on the vegetation–impervious surface–soil model

Abstract This study explores a spatiotemporal comparative analysis of urban agglomeration, comparing the Greater Toronto and Hamilton Area (GTHA) of Canada and the city of Tianjin in China. The vegetation–impervious surface–soil (V–I–S) model is used to quantify the ecological composition of urban/peri-urban environments with multitemporal Landsat images (3 stages, 18 scenes) and LULC data from 1985 to 2005. The support vector machine algorithm and several knowledge-based methods are applied to get the V–I–S component fractions at high accuracies. The statistical results show that the urban expansion in the GTHA occurred mainly between 1985 and 1999, and only two districts revealed increasing trends for impervious surfaces for the period from 1999 to 2005. In contrast, Tianjin has been experiencing rapid urban sprawl at all stages and this has been accelerating since 1999. The urban growth patterns in the GTHA evolved from a monocentric and dispersed pattern to a polycentric and aggregated pattern, while in Tianjin it changed from monocentric to polycentric. Central Tianjin has become more centralized, while most other municipal areas have developed dispersed patterns. The GTHA also has a higher level of greenery and a more balanced ecological environment than Tianjin. These differences in the two areas may play an important role in urban planning and decision-making in developing countries.


Spatiotemporal analysis of urban environment based on the vegetation-impervious surface-soil model 1 Introduction
Nowadays, cities across the world are spreading into their surrounding landscapes, sucking food, energy, water, and resources from the natural environment and profoundly changing global ecosystems. Urban ecosystems are among the most dramatic manifestations of anthropological impacts on the environment. 1,2 Recent research has also shown that urban sprawl and increased land use/land cover change in urban areas have considerable impacts on climate, hydrological, and biogeochemical cycles at various scales. 3,4 Therefore, it is important to understand clearly the different urban land surface characteristics, their spatial and temporal dynamics, and the corresponding responses to regional and global environment change.
The advantages of remote sensing-such as relatively low cost, large-area coverage, and repetitive observations-and also the growing interest in urban ecosystems have promoted the application of remote sensing in monitoring the urban environment. Researchers have developed a variety of approaches to distinguish surface cover types and numerous environmental variables through using remote sensing data in urban areas. Ridd 2 presented a concept model to parameterize the biophysical composition of urban environments. This model, namely the vegetation-impervious surface-soil (V-I-S) model, may serve as a foundation for characterizing urban/near-urban environments and has been applied in various studies. Chen 5 used Landsat TM and SPOT-P data to examine the effect of V-I-S cover types as input to hydrological models for runoff prediction for storms. Gluch et al. 6 utilized advanced thermal land applications sensor data to determine urban heat island effects using V-I-S cover types. Rashed et al. 7 employed four spectral bands of Indian IRS-1C imagery to study the anatomy of Greater Cairo in terms of endmember fractions (vegetation, impervious surface, soil, and shade) through linear spectral mixture analysis. Artificial neural network algorithms, [8][9][10] such as the multilayer perception feed-forward network and the back-propagation network, were applied for estimating end-member fractions within each pixel. Based on the V-I-S fractions, the classification and regression tree (CART) algorithm has also been widely used to obtain the impervious surface end-members. [11][12][13] More recently, some researchers have attempted a modern machine learning method-the support vector machine (SVM)-to map urban characteristics. [14][15][16][17] Meanwhile, studies of quantitative landscape ecology have developed landscape spatial metrics or indices to quantify the spatial contiguity, shape, and patch distribution of classified data.
A major problem is that, at the moment, it is difficult to extract accurately urban components based on remotely sensed data. Each image classifier has its own advantages and disadvantages. It is difficult to get a higher accuracy by relying on one classifier. Combining the strengths of several different methods-particularly through the synthesis of an advanced classifier and by using knowledge obtained from auxiliary data or empirical models-would contribute to producing a higher accuracy. In this paper, a machine learning methodology of global optimization -the SVM algorithm-was chosen as the classifier because of its improved performance compared to other approaches. This improved performance has been demonstrated in previous studies. [15][16][17] In addition, two kinds of knowledge were used to help to extract urban components: geospatial auxiliary data, such as urban vector data, and spectral index models, such as the modified normalized difference water index (MNDWI) 18 and normalized difference bareness index (NDBaI). 19 Although the V-I-S segmentation method of land cover is superior to models or classifications of land use for energy and moisture flux analysis of ecosystems, 20 there is still some spectral confusion because the same object can have different spectra or different objects can have the same spectrum. In this paper, four end-members including a high-albedo object, a low-albedo object, and easily distinguishable vegetation and soil were selected for the classification and fraction computation. After this had been done, empirical models and geospatial data were used to map this to V-I-S components.
Remotely sensed information has been widely applied in urban environment studies. Multitemporal remote sensing data have been used to map urban areas and reveal urban spatial and temporal changes. [20][21][22][23][24][25][26][27][28][29][30][31][32][23][24][25][26] demonstrates that the spatial metrics can be useful tools for monitoring a particular urban landscape by providing various measures of patch shapes and patterns. Most studies focus on case studies in developing countries [21][22][23]28,29 or in developed countries. 24,30 Moreover, comparison analyses between different countries are mostly based on coarse resolution imagery. 27,31 Due to difficulties in data acquisition and data processing, 6 not much research has been conducted into comparison analysis at medium resolutions. In 2011, the Comparative Study on Global Environmental Change Using Remote Sensing Technology project started. This program selects areas from China, Canada, Brazil, and Australia that are sensitive to global environmental change as test sites. Comparative studies of typical global change phenomena are performed between the northern and southern hemispheres, and also between the eastern and western hemispheres so as to develop the methodology of global environmental change remote sensing and to provide scientific data and decision support to deal with global environmental change. As a preliminary study of this project, we selected two urban areas, one in Canada and one in China, to explore the differences in urban development and the ecological implications of these differences.

Study Area
The Greater Toronto and Hamilton Area (GTHA) is located at the western end of Lake Ontario in Southern Ontario, Canada [ Fig. 1(a)]. Extending from 43°17′ N to 44°2′ N and 78°6′ W to 80°5′ W, it is the largest metropolitan area in Canada and is ranked in the top 50 most densely populated areas in the world. 33 The GTHA covers the city of Toronto as well as the surrounding regional municipalities of Durham, Halton, Peel, York, and Hamilton. The city of Tianjin is located between 38.57°N and 40.25°N and between 116.72°E and 118.32°E, covering ∼11;760.26 km 2 [ Fig. 1(b)]. Tianjin is one of the four directly governed municipalities of China and comprises the fifth largest area of urban land in China. In terms of urban population, Tianjin is the sixth largest city in China. Since 1980, the urban population of Tianjin has grown by 50% and has now reached 10.29 million. 34

Datasets and Preprocessing
Owing to their unique features of a long-time span, low cost, and relatively high temporal and spatial resolution, Landsat datasets were selected as our primary datasets in this study. The Landsat ETM+ and TM data were acquired from the global land cover facility at the University of Maryland (ftp://ftp.glcf.umd.edu/glcf/Landsat/), the U.S. Geological Survey website (http://eros.usgs.gov), and the remote sensing data sharing platform at the Center for Earth Observation and Digital Earth, Chinese Academy of Sciences (http://ids.ceode.ac.cn/). Four Landsat scenes (path/row of 18/29, 18/30, 17/29, and 17/30) were downloaded to cover the GTHA, while two scenes were downloaded (122/32 and 122/33) for Tianjin. Three stages of Landsat imagery from 1985 to 2005 were used to investigate the spatiotemporal dynamics of both study sites. To reduce classification uncertainty, only images that had limited cloud cover and that fell in the period of August to September (except for one 1987 scene of Tianjin) were selected. Detail descriptions of the datasets used are shown in Table 1. For the GTHA, three land use and land cover (LULC) maps (for 1985, 2001, and 2005) were provided by the Canadian Urban Land Use Survey. All Landsat scenes were processed to level 1T using the Universal Transverse Mercator map projection (UTM, Zone 17N for GTHA, UTM, Zone 50N for Tianjin), WGS84 geodetic datum, and north-up image orientation. All images and LULC maps were georegistered with RMS errors of <1 pixel. The LULC data were resampled to 10 m and the Landsat data were resampled to 30 m.

Methodology
In this study, a machine learning methodology of global optimization, the SVM algorithm, was chosen as the classifier because of its improved performance, as demonstrated in previous studies, compared with other approaches. 15,17 Moreover, two kinds of information were used to help to extract urban components: geospatial auxiliary data, such as urban vector data, and spectral index models. [18][19][20][21][22][23][24][25][26][27][28][29]35 Statistical analysis of fraction images and spatial metric analysis were then carried out in municipal districts and metropolitan areas in both study areas. The flow chart in Fig. 2 shows the processing chain used in this study.

Support Vector Machine
The SVM originates from a nonparametric machine learning methodology that follows basic ideas of supervised classification to map data into a high-dimensional space in order to identify the hyperplane of each class. [36][37][38] Due to its advanced generalization ability in targeting classification rules, it facilitates comparatively accurate computation and can carry out image classification processing quickly. 39 In this paper, about 800 samples for each stage image were manually selected for each endmember and were used for training. These included high-albedo objects, low-albedo objects, vegetation, and soil from the original images. Meanwhile, another 800 samples were selected and used for classification evaluation. The samples were selected from the most extreme areas of the feature space with the help of Envi4.7 software models (i.e., minimum noise fraction, pixel purity index, and n-dimensional visualization) and also using LULC category information. More particularly, the vegetation samples were taken from the areas of dense vegetation cover, the impervious ones from commercial/industrial areas or major road arteries, the soil from fallow agricultural land, and the water from areas of clear, deep water.
SVM uses training data to generate a model for estimating the targeted values. Training the SVM with a radial basis function requires two parameters: the penalty parameter controls the trade-off between maximizing the margin and minimizing the training error, while the gamma parameter describes the kernel width. In this study, the SVM library LIBSVM developed by Chang and Lin 40 is employed. The influence of both parameters on the producer's accuracy (PA) is analyzed in the process. In Fig. 3, the classification accuracy improves sharply at a gamma parameter of 0.1 and changes slightly in the range [0.1, 0.6]. The SVM performs better with a penalty parameter in the range [100, 300] and its accuracy decreases at larger values. To achieve a better accuracy and a faster processing time, in this study, the gamma parameter and penalty parameter were set to 0.143 and 100, respectively.

V-I-S and Knowledge-Based Mapping Models
The V-I-S model [Eq. (1)] is a conceptual framework that groups the great variety of urban land covers into three general categories-vegetation, impervious surface, and soil-plus water. 1 ρ where ρ i , ρ imp;i , ρ veg;i , ρ soil;i , and ρ water;i are the reflectance for land type, impervious surface, vegetation, soil, and water, respectively, for band i; f imp , f veg , f soil , and f water are the corresponding fractions of associated components; and e i is the residual error.
Previous studies 35,41 demonstrated that impervious surfaces were located on or near the line connecting the low-albedo and high-albedo end-members in the feature space. Therefore, the percentage covered by impervious surfaces can be estimated as a linear combination of lowalbedo and high-albedo surfaces. Then, Eq. (1) can be described as where ρ i ,ρ low;i , ρ high;i , ρ veg;i , and ρ soil;i are the reflectance for land type, low-albedo, high-albedo surface, vegetation, and soil, respectively for band i; f low , f high , f veg , and f soil are the corresponding fractions of associated components; and e i is the residual error. Equation (2) shows that spectral reflectance of an urban surface is composed of four endmembers (high-albedo objects, low-albedo objects, vegetation, and soil). The mapping model between V-I-S and the end-members was established based on Eqs. (1) and (2) and on knowledge-based empirical values of the indexes shown in Fig. 2. Generally, low-albedo materials comprise water, dark impervious surfaces (i.e., asphalt road), and shadow. Water could be easily extracted according to the low-albedo classification result and the empirical value of MNDWI (MNDWI > 0.1) derived using Eq. (3). 18 Other information such as geospatial vector data (for example, urban administrative boundaries and road network data) is very useful for separating different components that have similar spectra. Shadows in urban area make up only small proportions of the whole image and so can be ignored. Areas in shadow outside urban areas generally consist of the shaded slopes of hills, which can be separated into vegetation or soil by using the NDVI [Eq. (4)]-threshold for soil: NDVI < 0.2; otherwise, vegetation. High-albedo materials such as dry soil and sand tend to be confused with bright impervious surfaces. Highalbedo impervious surfaces can be extracted by removing dry sands or dry soils with NDBaI [Eq. (5)] (dry soil: NDBaI > 0). Following the steps in Fig. 2, we can get four fraction images of V-I-S and one V-I-S classification image.

Spatial and Temporal Analysis Models
In this study, we used two indicators or models to measure the spatial-temporal dynamics of the urban environment: rate of change of impervious surface and the aggregation index (AI).

Rate of change of impervious surface
This is calculated using Eq. (6) and compares the dynamic differences of urban impervious surfaces. Fig. 3 Three-dimensional graph of gamma parameter (x ), penalty parameter (y ), and producer's accuracy (z).
where C i is the rate of change of impervious surface during the study period, W ei and W bi are impervious surface area or percentage at the end and beginning of the study period, respectively, and t is the length of the study period measured in years. The results of this calculation give the annual rate of change of the impervious surface.

Aggregation index
The study also provided quantitative measures of urban growth pattern besides extent and change estimates in the two urban areas during the two stages. In the field of landscape ecology, a number of metrics have been developed for assessing urban dispersion. These include the landscape shape index, AI, split index, and so on. In this research, the AI [Eq. (7)] was selected to measure class aggregation or clumping, which increases as the patch type becomes more aggregated. The index is easily calculated using ArcGIS or Fragstats software.
where g ii is the number of like adjacencies between pixels of patch type i based on the singlecount method, max -g ii is the maximum number of like adjacencies between pixels of patch type i based on the single-count method, and P i is the proportion of landscape composed of patch type i. In this study, the V-I-S outputs were aggregated into six classes according to impervious surface percent (ISP), that is, high albedo (ISP [0.8,1]), middle albedo (ISP [0.5,0.8]), low albedo (ISP [0,0.5]), vegetation, and soil. Then we explored the zonal spatial patterns of the six classes in the study area (as shown in Fig. 1) using a spatial metric analysis software (FRAGSTATS 3.3).

SVM Classification
Four classes-high albedo, low albedo, vegetation, and soil-were derived from Landsat images at each stage. Accuracy assessment was conducted only for the GTHA as only GTHA has matched evaluation data. For accuracy assessment, the confusion matrix approach was applied to extract four accuracy indicators: overall accuracy (OA), user's accuracy, PA, and the overall kappa coefficient (OK). The confusion matrix is shown as Table 2.
From the table, it can be seen that, in this study, this type of segmentation method produced good results with a high OA (96.69%) and high OK (95.58%). The relatively lower separability between high albedo and soil resulted in a lower PA (86.96%) compared to other end-members.

VIS Fraction Components
Using the V-I-S model, four fraction images-vegetation, impervious surface, soil, and waterwere developed for Tianjin [ Fig. 4(a)] and the GTHA [ Fig. 4(b)] in three stages. Three statistical indicators comparing the estimated V-I-S fraction and the reference data [root-mean-square error (RMSE), the mean absolute error (MAE), and the coefficient of determination (R 2 )] were used to evaluate the accuracy of the V-I-S estimation. Detailed formulae and descriptions of the assessment indicators can be found in Ref. 42. One hundred pixels were able to be selected using a random function from LULC (Fig. 5) data (only data from 1985 and 2005 were selected because of considerations of synchronicity), and based on the selected pixels, one hundred 90 m × 90 m samples for statistical analysis were produced (Fig. 6). Due to the different categories in the V-I-S land surface cover and the LULC data, the classification results could not be directly compared with LULC data. For this reason, knowledge-based indexes were used to build a set of mapping rules for comparison (Table 3). In Table 3, the numbers in brackets [(1) to (9) (5); case (10): for transitional urban recreational areas and other LULC, discard the sample and then reselect. The urban component fractions equal the sum of the corresponding weights of each category weight divided by the number of effective pixels (9 pixels × 9 pixels − discarded pixels). The error statistics of the total fractional covers were calculated by comparing the V-I-S mapping model with the LULC data for the GTHA (Table 4).  The low RMSE and MAE and high R 2 of all the fraction components also indicated the high accuracy of the approach. Among popularly used subpixel classifiers, the linar spectral mixture analysis-based method tends to overestimate the impervious surface fraction in nonurban areas and to underestimate it in urbanized areas. 9 CART-based methods are more sensitive to data noise. 13 Impervious surface is often confused with dry soil and water. In this case, it is difficult to accurately segment impervious surface from other objects directly. In our research, we selected the four end-members according to spectral characteristics and took advantage of the comparatively accurate and fast processing capability of the SVM algorithm and knowledge-based empirical models to obtain a higher level of accuracy.

Comparison of Spatiotemporal Patterns in GTHA and Tianjin
Using the V-I-S fraction components, the spatial and temporal dynamics of each component can be easily depicted. The V-I-S components in municipal districts and metropolitan areas were extracted from the fraction images. The two graphs-one for the GTHA [ Fig. 7(a)] and one for Tianjin [ Fig. 7(b)]-reveal the patterns and processes of urban development in both areas. In the past 20 years, urban areas have experienced rapid growth. However, there is a great difference between the situations in the GTHA and Tianjin. For the GTHA, it mainly increased from 1985 to 1999. After 1999, only two municipal areas (Peel and Halton) and four cities (Brampton and Mississauga in Peel, Vaughan in York, and Burlington in Halton) have positive urban growth trends. In 1985, only Toronto had a high level of urbanization, while the other areas had relatively low urbanization levels. Thus, urban development in the GTHA could be regarded as having a monocentric pattern. After 1999, satellite cities or towns adjacent to Toronto, including Brampton and Mississauga in Peel and Vaughan in York, reached high levels of urbanization and became new urban centers. At the same time, other satellite cities or towns farther away from Toronto became subcenters. This could be called a multicentric urban pattern. Tianjin has experienced continuous rapid urban sprawl since 1985. Since 1999, the rate of urban sprawl has accelerated at the cost of soil and vegetation loss. In 1985, only the central district had a high level of urbanization, but in recent decades, the coastal area has developed into a new center (Binhai New Area). In the ternary diagram (Fig. 8) of the urban environment, all sites in the  GTHA are near to the vegetation and impervious surface axes, while the ones in Tianjin are near to the soil (plus water) and impervious axes, which demonstrates the great differences between the two areas. In the 1980s, soil characteristics are clear for Tianjin due to the image acquisition date (May); impervious surface fractions at most sites in the GTHA are higher than in Tianjin in the 1980s and in 1999 but are lower in 2005. In addition, the impervious surface fraction for Tianjin is increasing at an accelerating rate all the time, while for the GTHA a point of inflexion has been passed. Finally, the GTHA has a high level of urban greenery, while in Tianjin it is very low. Figure 9 shows the annual rate of change in the impervious surface area for the GTHA and Tianjin. Figure 9(a) indicates that from 1985 to 1999, all sites in the GTHA had positive change rates, whereas from 1999 to 2005, the rates of change at some sites had turned negative. It can also be seen that in the first stage, York municipality and its city Vaughan had a high rate of change due to its proximity to Toronto; in the second stage, the rate of change at all sites was far slower than in the first stage. From Fig. 9(b), it can be seen that in the first stage, the rate of change was ∼5%, whereas in the second stage, it had reached 25% in the Beichen district. Figure 10 illustrates the spatiotemporal dynamics of the urban growth patterns of the study areas. From Fig. 10(a), it can be seen that all districts in the GTHA became more disaggregated from 1985 to 1999 and more aggregated from 1999 to 2005. For Tianjin [ Fig. 10(b)], it can be seen that more patches became aggregated in the first stage and became more decentralized in urban areas in the second stage, except for central Tianjin.

Discussion
In this study, we used the MNDWI to remove the water from the low-albedo fraction image and used the NDBaI to eliminate dry soils from the high-albedo fraction image. The evaluation of the accuracy showed that the accuracy was high. However, we could not remove the shade from urban areas in the low-albedo fraction image, although this could correspond to impervious surface. Also, the NDBaI index could only eliminate very dry soils (primary barren soil). In the evaluation process, samples containing transitional and urban recreational areas were discarded. In fact, these areas could contain multiple categories. Therefore, we need to evaluate the validity of this method in future work, focusing on transitional urban areas. In addition, this assessment was carried out for the GTHA only due to the lack of synchronous high-resolution images for Tianjin.   The spatiotemporal dynamics of the GTHA and Tianjin revealed that there were great differences between the two urban areas. This could be related to factors such as population and policy. According to census data, the proportion of the population classified as urban in Tianjin was 53.70, 57. 29, and 75.11% in 1985, 1999, and 2005, respectively. In the GTHA, it exceeded 82% in 1985. Figure 11(a) illustrates population and population density in Tianjin, the GTHA, and their central cities (central Tianjin and Toronto). The most surprising feature of this graph is that the population density of central Tianjin far exceeds that of the other areas, although there is no obvious population density difference between Tianjin and the GTHA and other pairs of population and population density curves (GTHA and Tianjin; Toronto and central Tianjin) are almost parallel. The large increase in the urban population percentage in Tianjin and the high population density in central Tianjin may partly account for the rapid increase in impervious surfaces. However, population and population density increases in the GTHA and Toronto did not bring about a corresponding appreciable increase in impervious surfaces from 1999 to 2005, which implies that another factor played an important role in the urban development. Since the 1990s, experts in urban studies and environmentalists in Canada have raised concerns about the region's rapid, low-density pattern of land development, calling for better management of future land-use decisions. The provincial government has also introduced related legislation, such as the "Places to Grow Act 2005." In contrast, China experienced a real estate bubble in the 1990s, which resulted in a disorganized expansion of developed land. These different policy principles are another factor impacting on the dynamics of the urban ecological environment.

Conclusions
This study explores the application of the V-I-S method in large-area urban agglomeration mapping and ecological environment analysis. First, four end-members (high-albedo surface, lowalbedo surface, vegetation, and soil) were selected for classification using the SVM algorithm. The accuracy of the classification results reached 96.69%. Next, to map the SVM results to V-I-S component fractions, a knowledge-based method was applied to the results to establish a mapping model-the RMSE ranged from 3.12 to 6.43% and R 2 ranged from 0.94 to 0.99. Finally, a spatiotemporal analysis was conducted for municipal and metropolitan areas in Tianjin, China, and the GTHA, Canada. We found that both urban areas experienced rapid growth in the period from 1985 to 2005. For the GTHA, urban sprawl mainly occurred from 1985 to 1999. Urban development patterns underwent a change from monocentric (based on Toronto) to multicentric (based on Toronto, Brampton, and Vaughan). For Tianjin, all districts underwent rapid urban expansion, with the rate accelerating after 1999. The aggregation index results reveal that the GTHA and its municipal areas became more disaggregated from 1985 to 1999 and then more aggregated from 1999 to 2005, which implies that the Canada government's initiatives achieved their intended result of compact growth. Most urban sites in Tianjin showed the opposite trend. The results reveal different patterns of sprawl in the GTHA and Tianjin, which may reflect the demographic and economic influences on urban environments across the globe.
Qingni Huang received her PhD from Center for Earth Observation and Digital Earth, Chinese Academy Science in 2012 and is now working in National Satellite Meteorological Center. Her research interests include urban environment, retrieval and simulation of land surface temperature.
Xinwu Li is an associate professor in the Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences. He has published 30 papers about polarimetric and interferometric SAR remote sensing. His current research interest is global environment change especially in urban area based on remote sensing. Zhongchang Sun is an assistant professor of Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences. He has issued three science citation index papers about urban impervious based on remote sensing.
Ying Zhang is a research scientist at the Canada Centre for Remote Sensing, Natural Resources Canada. She has been working on research and technical development of optical remote sensing and atmospheric physics since 1987. Her current research interests include applications of highresolution remote sensing to urban mapping, disaster management, and sustainability issues related to natural resource development.