2.1.

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.³³ 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 $\sim \mathrm{11,760.26}{\mathrm{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.³⁴

Fig. 1

Location map of the study areas is (a) Greater Toronto and Hamilton Area (GTHA) and (b) Tianjin.

2.2.

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\text{pixel}$. The LULC data were resampled to 10 m and the Landsat data were resampled to 30 m.

Table 1

Characteristics of satellite data used for the study area.

3.1.

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–38 Due to its advanced generalization ability in targeting classification rules, it facilitates comparatively accurate computation and can carry out image classification processing quickly.³⁹

In this paper, about 800 samples for each stage image were manually selected for each end-member 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⁴⁰ 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.

Fig. 3

Three-dimensional graph of gamma parameter ($x$), penalty parameter ($y$), and producer’s accuracy ($z$).

3.2.

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.¹

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 low-albedo and high-albedo surfaces. Then, Eq. (1) can be described as

Equation (2) shows that spectral reflectance of an urban surface is composed of four end-members (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 ($\mathrm{MNDWI}>0.1$) derived using Eq. (3).¹⁸ 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: $\mathrm{NDVI}<0.2$; otherwise, vegetation. High-albedo materials such as dry soil and sand tend to be confused with bright impervious surfaces. High-albedo impervious surfaces can be extracted by removing dry sands or dry soils with NDBaI [Eq. (5)] (dry soil: $\mathrm{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.

3.3.

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).

4.1.

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.

Table 2

Confusion matrix of classification accuracies.

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.

4.2.

VIS Fraction Components

Using the V–I–S model, four fraction images—vegetation, impervious surface, soil, and water—were 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\mathrm{m}\times 90\mathrm{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)] represent different special circumstances as determined by the LULC data, the knowledge-based indexes (NDVI and MNDWI), and empirical threshold values:^18,43 case (1): $\mathrm{NDVI}>0.5$, $\text{weight}=1$; case (2): $\mathrm{NDVI}<0.2$, $\text{weight}=0$; case (3): $0.2\le \mathrm{NDVI}\le 0.5$, discard the pixel; case (4): $\mathrm{NDVI}>0.5$, $\text{weight}=0$; case (5): $\mathrm{NDVI}<0.2$, $\text{weight}=1$; case (6): $\mathrm{MNDWI}>0.1$, $\text{weight}=1$; case (7): $\mathrm{MNDWI}<0.1$, $\text{weight}=0$; case (8): $\mathrm{MNDWI}<0.1$ and (2), $\text{weight}=0$; case (9): $\mathrm{MNDWI}<0.1$ and (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\text{pixels}\times 9\text{pixels}-\text{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).

Fig. 4

Vegetation, impervious surface, soil, and water fraction images: (a) Tianjin and (b) GTHA.

Fig. 5

Land use and land cover data from the Canadian Urban Land Use Survey.

Fig. 6

Accuracy assessment flow chart from the vegetation–impervious surface–soil (V–I–S) mapping model.

Table 3

Mapping rules for mapping weights of land use and land cover (LULC) contribution to vegetation–impervious surface–soil (V–I–S) fraction (1 means 100%; 0 means 0%).

Table 4

Accuracy assessment of fractional estimation results for V–I–S components.

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.⁹ CART-based methods are more sensitive to data noise.¹³ 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.

4.3.

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.

Fig. 7

The temporal–spatial dynamics of VIS components: (a) GTHA and (b) Tianjin.

Fig. 8

A dynamic ternary V–I–S diagram for municipal districts and metropolitan areas in the GTHA and Tianjin areas.

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 $\sim 5\%$, whereas in the second stage, it had reached 25% in the Beichen district.

Fig. 9

Annual rate of change rate of impervious surfaces: (a) Tianjin and (b) GTHA.

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.

Fig. 10

The spatiotemporal dynamics of urban dispersion in municipal districts and metropolitan areas: (a) GTHA and (b) Tianjin.

4.4.

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.

Fig. 11

Population time series.