Urban growth sustainability of Islamabad, Pakistan, over the last 3 decades: a perspective based on object-based backdating change detection

Urban growth copes with problems in sustainable development. In developing countries, particularly, sustainable development of urban growth copes with severe challenges with respect to sluggish economic and social growth, population boom, environmental deterioration, unemployment, slums and so on. Time series of remote sensing data provide critical support on sustainability assessment. However, the urban spatial extend cannot be accurately extracted from land cover data. Targeting the urban growth and its sustainability in Islamabad, the capital of Pakistan, this study extracts urban area from four periods of Landsat images between 1990 and 2018 using an innovative object-based backdating change detection method and two criteria for extracting urban land from impervious surface. We prove that impervious surface cover and urban area increased 273.10% and 426.21%, respectively, over the last 3 decades. We identify five factors playing important role in urban growth: population, transportation systems, master planning, industrial and real estate development, and neighbor urban effect. In this study, we assess the socio-economic sustainability associated with slum growth and census data, and the environmental sustainability in relation to the variations of normalized difference vegetation index (NDVI) in forest areas. We found that slums increased with the corresponding growth of urban area and population, reflecting sluggish economic increase in Islamabad. We found that the area of woodland increased 9.29%, but its NDVI decreased from 0.668 to 0.551, implying a deteriorative trend of environmental condition.


Introduction
Urban growth is an important trend due to the needs of immigrants who move from other parts of the country. Until 2018, more than 55% of the world population lived in cities. This proportion is expected to reach 68% by 2050 (United_Nations 2018). Rapid urban growth normally exceeds the capacity for local governments to deliver services and infrastructure, which increases urban poverty and slums, especially in developing countries (Duque et al. 2015). Urban growth further leads to serious problems in cities, including urban flooding, pollution, epidemic diseases, social inequality, and political instability (Netzband et al. 2007). Thus, the study of urban growth sustainability is of great significance for decision making about urban planning in relation to construction programs, poverty reduction, disaster monitoring and prevention, and improvement of the environmental quality (Bhatta 2010). Current urban sustainability research is mostly based on field socio-economic surveys and demographic statistics, which is not only time consuming and labor intensive, but also makes it difficult to grasp spatial heterogeneity within cities.
Remote sensing technologies provide necessary data support in the study of urban growth sustainability. Medium-resolution remote sensing images, such as Landsat Thematic Mapper images with a spatial resolution of 30 m, allow 16-day interval coverages for monitoring global land cover changes since the 1980s (Woodcock et al. 2008;Zhu et al. 2019). Previous studies generally used pixel-based image analysis to extract the urban area (Hassan et al. 2016;Butt et al. 2012Butt et al. , 2016Mannan et al. 2018). Most studies are less concerned about the inner structure of the city because of classification error associated with the so-called ''salt and pepper effect'' owing to the method (Blaschke 2010;Blaschke et al. 2000Blaschke et al. , 2014Durieux et al. 2008). Further, the urban spatial extend cannot be of course extracted from land cover data (Li et al. 2016). Many studies obtain impervious surfaces from remote sensing images and assume that they correspond to urban areas (Weng 2011;Lu et al. 2010), however, impervious surfaces not only exist within cities but also in various settlements in suburban and rural areas.
In this study, we selected Islamabad City, the capital of Pakistan, for assessing urban growth sustainability. We improve the image analysis approach with an innovative object-based backdating change detection (OBBCD) method to extract impervious surface. Further, we introduce two spatial criteria to distinguish the urban area from impervious surface. We explore the socio-economic sustainability of Islamabad through the analysis of the development and evolution of slums. We use normalized difference vegetation index (NDVI) as a measure of forest ecosystem health to evaluate environmental sustainability.
Following the introduction, ''Literature review'' section reviews literatures about urban growth assessment based on remote sensing image analysis. ''Study area and data preparations'' section presents the conditions of Islamabad city and data preparation. ''Method'' section addresses the OBBCD approach, spatial criteria for extracting urban areas from impervious surface and sustainability assessment based on remote sensing technology. ''Results'' section evaluates the land cover classification, presented transition to impervious surface from other land cover types, and proposed growing patterns of urban area. ''Discussion'' section discusses factors affecting urban growth in Islamabad and socio-economic and environmental sustainability. Finally, ''Conclusions'' section summarizes the research.

Literature review
Remote sensing is an advantageous earth observing means for monitoring and assessing urban sustainability (Rochon et al. 2004;Trinder 2017). The widely used earth observing data are moderate-resolution Landsat images, particularly Landsat 1/2/3 MultiSpectral Scanner (MSS), Landsat 4/5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper plus (ETM?) and Landsat 8 Operational Land Imager (OLI) often at no cost (Hansen and Loveland 2012;Zhu 2017;Zhu et al. 2019). Other commonly used moderate-resolution imagery include SPOT (Systéme Probatoire l'Observation de la Terre), ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) and air photos. Although high-resolution remote sensing technology with spatial resolution less than 10 m has advanced rapidly in the last 20 years, moderate-resolution images dominated monitoring and assessing methods are remain irreplaceable because of its large periodic coverage for nearly half a century. The first step for monitoring and assessing urban sustainability is extracting land cover and its changes. Most studies used Maximum Likelihood Classifier (MLC) combined with post-classification change detection. With this method and moderateresolution data, Dewan and Yamaguchi (2009) Hepcan et al. (2013) studied the urban growth of Izmir, Turkey; Moghadam and Helbich (2013) established the relationship between the expansion of urban built-up area and the retreat of open land and cropland in Mumbai, India, from 1973Boori et al. (2015) found the land cover change and corresponding population growth and decline in Samara, Russia, from 1972Russia, from to 2009Hassan et al. (2016) showed that cropland, built-up areas and water bodies increased, and forest and bare land decreased in Islamabad from 1992 to 2012; Pourebrahim et al. (2015) revealed land cover changes in Kuala Langat district, Malaysia, from 1988 to 2010 and predicted urban growth pattern in 2025; Morshed et al. (2017) indicated that the sprawl of Dhaka City, Bangladesh, were relevant to the shrinking of natural vegetation, farmland and water body from 1989 to 2014; Fenta et al. (2017) showed the relationship between the growth rate of built-up area and land cover type conversion in Mekelle, Ethiopia, from 1984Sandamali et al. (2018) studied the urban land cover change in Kuala Langat district, Malaysia, from 1980 to 2010; Enoguanbhor et al. (2019) found that the urban growth in Abuja, Nigeria, was at the cost of the reduction of vegetation area from 1987 to 2017. Besides, few studies tried new methods. Xi et al. (2018) adopted random forest classifier and postclassification change detection for extracting the urban area in Xuzhou, China, from 1985to 2015. Dereli (2018 revealed urban land use change in Istanbul, Turkey, from 2003 to 2016, by using multi-layer perceptron artificial neural network. Aina et al. (2019) compared spectral angle mapping, spectral mixture analysis and band ratioing for land-cover classification in Riyadh Saudi Arabia, from 1972 to 2014.
Nevertheless, the above studies are generally based on pixel-based spectral analysis. Owing to the effect of ''different objects performing with similar spectral features and similar objects showing diverse spectral features'' in images associated with the pixel-based approach, as well as the so-called ''pepper and salt'' effect, there are inevitably serious errors in the results of land cover classification and change detection (Blaschke et al. 2000;Radoux and Bogaert 2017). In the last 2 decades, object-based image analysis (OBIA) has been developed. As the pre-requisite step, image segmentation divides an image into internally homogeneous, mutually disjoint segment objects (i.e., polygons), and assigns all pixels to separate segment objects. The sparkle noise which leads to the ''salt and pepper effect'' in pixel-based image analysis is thus smoothed into segment objects. A large number of studies proved that OBIA is superior to the pixel-based approach, owing to its efficient adoption of spectral, textural, geometric and contextual features of segment objects in classification (Blaschke 2010;Blaschke et al. 2014;Cheng and Han 2016;Chen et al. 2018a, b;Myint et al. 2011). Furthermore, a novel object-based change detection (OBCD) method has been developed by combining OBIA with post-classification change detection (Chen et al. 2012). Esbah et al. (2009) used OBCD method for analyzing land cover changes in Didim peninsula, Turkey. (Keshtkar et al. 2017) studied the conversion of urban land cover types from 1990 to 2010 in Thuringia, German, using an objectbased support vector machine classifier combined with change detection method Shen et al. (2017). explored urban growth from 1979 to 2013 in Chongming Island in Shanghai, China. On the other hand, a lot of work remains to improve the urban land cover change approach by making sufficient use of medium-resolution image for a long time in OBCD method.

Study area
Islamabad is located at the southeast piedmont plain of Margalla Hills, with altitude ranging between 457 and 1240 m above sea level (Fig. 1). The region has a tropical/subtropical mountain monsoon climate with mean annual temperature of 20.9°C and annual precipitation that reaches 1323.4 mm in plain area. The dominant natural vegetation is tropical evergreen broad-leaf forests at an altitude of less than 900 m and subtropical evergreen coniferous forest and deciduous broad-leaf forest above this elevation (Harry and Seth 1965). With the intensification of human activities, most natural forests have moved and succeeded by secondary shrubbery, grass, and crops. Cultivation is divided in two seasons. Winter crops are sown during October and December and harvested during March and April, while summer crops are grown from February to October.
Islamabad City was built when the Pakistan government moved the capital from Karachi in the 1960s. Since then, the population has multiplied. In 1961, there were only 118 thousand people in Islamabad. The population raised to 237.5 thousand in 1972, 520.2 thousand in 1992, 1.22 million in 2009 and more than 2 million in 2017. People from all over the country flocked to the capital, hoping to take advantage of its development to improve their employment, living conditions and education.
The administrative region of Islamabad is divided into five zones (Fig. 1). The master plan was limited in both Zone-I and -II (Butt et al. 2012) and divided the plan construction area into sectors with a size of about 2 km 9 2 km. They are given names from north to south from A to I and is further numbered from east to west in numerical order from G1 to G18. Some sectors are developed well, while others are under-developed. Zone-III is characterized by mountain forests and piedmonts. Most area in Zone-III refers to the Margalla Hills national park and is endowed with ecological protection. Zone-IV and the north Zone-V are mainly characterized by agricultural landscape, while the south Zone-V close to Rawalpindi City is characterized by private houses, social and industrial facilities.
Zone one and two in Islamabad were built strictly according to the master plan. The main roads in urban area interweave vertically and horizontally, dividing the whole urban area into dozens of districts of roughly equal size. They are functionally divided into administrative, public utilities, central business, residential, colleges and universities, industrial, foreign embassies, green districts and so on. Buildings in the administrative district include the presidential palace, the parliament building, the offices of government ministries, and foreign embassies and consulates. Commerce, transportation, communication and municipal construction are all well developed. The residential area is crisscrossed by streets that divide it into 76 square blocks, each with elementary schools, middle schools, recreational centers, shops, restaurants, mosques and other living facilities. Colleges and universities include Quaid-e-Azam University, Allama Iqbal Open University and Jinnah University. Industry is limited to small and medium-sized light industries such as food, furniture, construction and printing, which are closely related to the daily life of the residents. The agricultural areas, located outside the planned urban areas, mainly develop poultry and vegetable production.
Local reports showed that the environmental protection program has been carried out in the Margalla Hills national park. Under the guidance of the Capital Development Authority (CDA), which was founded in 1960 by the Pakistan government, hundreds of thousands of trees has been planted on Margalla Hills national park in the last 2 decades (https://www.cda. gov.pk). Another non-governmental organization, Margalla Hills Society has planted more than twenty thousand trees in the national park, including pine, Bauhinia variegata, and bottle brush (Callistemon).

Data preparation
We collected four Landsat TM/ OLI images with path 150 and row 37 from the United States Geological Survey (https://glovis.usgs.gov/) ( Table 1). These images were acquired in 1990, 1998, 2009, and 2018. Hereafter, we named them as IM1990, IM1998, IM2009, and IM2018, respectively. The four images were captured during February and March because the two months represent the most vigorous season for winter/spring crops in the study area. This is conducive to distinguish impervious surfaces from other land cover types. These images were characterized by zero cloud coverage in the study area. Then, we used the software package ENVI FLAASH (Excelis Visual Information Solutions Inc, Boulder, Colorado) to make radiometric calibrations and atmospheric corrections.
In data collection, we derived basic geographic data, such as place names, streams, lakes and reservoirs, roads, railways, airports, and slum areas by means of Google Earth. Field investigation was finished by one of the authors.
We applied a taxonomy system of land cover based on the Vegetation-Impervious surface-Soil (VIS) model (Ridd 1995;Weng and Lu 2009;Gluch and Ridd 2010). We divided ''vegetation'' in this model into woodland, cropland, shrub and grassland; renamed ''soil'' as bare land; and appended a ubiquitous type of water body. In practice, bare land could be fallow or temporary construction land. Thus, the taxonomy system consists of six types: cropland, woodland, shrub and grassland, bare land, water body, and impervious surface. Furthermore, we combined with Google Earth to complete the selection of sample data of various land cover types, including change and not-change parcels, which are respectively used for parameter optimization in image segmentation and classification, and accuracy assessment of classification results.
To improve the classification accuracy in impervious surface, we also used the Defense Meteorological Satellite Program/Operational Linescan System (DMSP/OLS) and the Suomi National Polar-orbiting Partnership/Visible Infrared Imaging Radiometer (NPP/VIIRS) night-light images in these four dates with respect to these Landsat images (Table 1). We performed reprojection, resampling, and relative radiometric correction on the night light data, according to the spatial resolution and coordinate system of Landsat TM/OLI images.

Object-based backdating change detection (OBBCD)
Many change detection methods have been developed for extracting land cover types and their changes from remote sensing images. The common change detection methods can be categorized into direct comparison, transformation, and post-classification approach (Hussain et al. 2013;Tewkesbury et al. 2015). Typical examples of direct comparison method use difference value, ratio, or regressive factors of a corresponding spectral channel in time series of multispectral images. Image difference or ratio conveniently supplies change/non-change information and while they are advantageous on operational simplicity, the deficiency lies in the lack of information about the land cover types and ''from-to'' change information among them. In addition, it is often problematic to differentiate changes owing to either land cover type or natural phenology. The transformation method is basically consistent with the direct comparison method. The prominent difference between them is that the transformation method requires multispectral operations of the time series of images, such as vegetation indices, principal component transformation, and change vector analysis. The disadvantage of the transformation method is similar to direct comparison: a lack of type change information.
The post-classification method is firstly based on image classification of two periods and then carries out overlaying operation on the classification data of different periods. It has the advantage of obtaining ''from-to'' change results, but the disadvantage is that the error in classification could propagate and magnify with increasing processing steps (Burrough and McDonnell 1998).
In recent years, a backdating method has been proposed (Linke et al. 2009). The method derives a pixel-based classification result from a recent image at first, detects pixel-by-pixel changes between recent and previous images afterwards, and then classifies the change pixels in previous images. The land cover types of non-change pixels are inherited from the recent image, while the ''from-to'' relationship of the change pixels can be obtained according to the results of two classifications.
Integration between OBIA and backdating change detection demonstrated high efficiency and high accuracy (Yu et al. 2016;Toure et al. 2018). Based on image segmentation and the segment object classification of a recent image, backdating analysis can distinguish change and non-change objects through change detection and threshold analysis. Each land cover type in recent image is assigned to the objects that did not change in previous image. For the change objects, their types in the previous image can be identified by another classification process. The ''from-to'' relationship of segment objects can be derived accordingly. Then, the OBIA results are GISready data with distinct boundaries of each land cover parcel (Blaschke 2010;Jalan 2012;Benz et al. 2004).
In this study, we first obtained segment objects from IM2018 under the support of the fractal network evolution algorithm (FNEA) (Baatz and Schäpe 2000) in the software eCognition Developer. We used an optimal combination of scale factor of 12, shape factor of 0.9, and compactness factor of 0.4 as defined in FNEA to achieve an optimal segmentation based on a trial-and-error process.
Second, we used a random forest algorithm (Breiman 2001) in segment classification. Before this operation, we selected the best 16 from a total of 70 candidate spectral and textural features associated with segment objects (see Table 2). To enhance the information about vegetation, building, water, and residential area, four spectral indices (NDVI, MNDWI, NDBI, and RRI) associated with the brightness, greenness and wetness components from Tasseled Cap Transformation (TCT) were considered in these candidate features. The calculation about TCT brightness, greenness, and wetness follows (Liu et al. 2015). The calculations about grey-level cooccurrence matrix (GLCM)-based textural features can be referred in the eCognition reference book (Trimble 2014). After optimization, 70 decision trees with unlimited depth and ten-fold cross-validation were adopted in this process. Thus, we obtained the classification result. Later, we called it LD2018.
Third, we used a fuzzy c-means clustering algorithm (Ghosh et al. 2011) to differentiate the change objects and no-change objects between IM2018 and IM1990, between IM2018 and IM1998, and between IM2018 and IM2009. In this process, we obtained the difference segment objects from the recent and previous images, from which the change and non-change objects were distinguished. The memberships of the segment objects belonging to change and nonchange clustering centers were measured using minimal Euclidean Distance. Further, we adopted the support vector machine algorithm (Hussain et al. 2013;Otukei and Blaschke 2010;Petropoulos et al. 2012) for classification of change objects to identify their categories at the previous period. When running this algorithm, a radial basis function was adopted as the kernel function, and the penalty coefficient and gamma coefficient were set as 2 and 0.08, respectively. Finally, we constructed the ''from-to'' schemes of land cover transition based on change objects and obtained the classification results in 1990, 1998, and 2009, namely LD1990, LD1998 and LD2009.
Extracting the urban area from remote sensing images Here, we regard cities and towns as urban area. It can be extracted from impervious surfaces. Impervious surface is generally regarded as a special land cover type related to buildings, roads, sidewalks, parking lots, and other urban objects (Weng 2011). It can be derived from land cover data. But impervious surface does not simply correspond to urban area because it could involve all sizes of settlements, including villages in countryside. It is therefore necessary to establish a set of criteria for extracting urban area from impervious surfaces and removing rural parts.
Considering urban areas as a subset of impervious surfaces, we proposed two criteria for extracting physical urban area from impervious surface data. The first criterion may be called the minimum area criterion: an impervious surface patch is regarded as urban area if its area is more than 0.1 km 2 . The size of an impervious surface patch reflects agglomeration degree of land used for urban residential use, commerce, or industries and transportation. Over small unit areas are not enough to accommodate people in diverse jobs, thus not leading to development of public service systems in an urban area. This study shows that Textural features of segment objects GLCM mean in all directions See eCognition reference book (Trimble 2014) GLCM correlation at the direction of 45°G LCM dissimilarity at the direction of 0°G LCM dissimilarity at all directions the proposed area threshold is reasonable in Islamabad.
The second criterion can be named as the proximity criterion. The impervious surface patch that extends outward of the urban area less than 100-m is concerned as a part of urban area. This distance was used to differentiate peri-urban from rural areas (Dupuy et al. 2012). This distance limit usually meets the requirements of urban transportation accessibility.

Sustainability assessment
The study of sustainable development mainly involves social, economic, and environmental aspects. The available method for sustainability assessment is based on statistical socio-economic, environmental and census data from ground investigations. In developing countries, however, socio-economic statistical and census data are often fragmentary, impeding the assessment of socio-economic sustainability (Stoler et al. 2012). Meanwhile, remote sensing observation provide indirect evidences for assessing urban sustainability. The socio-economic sustainability in urban area can be indirectly reflected by the difference in urban growth patterns. When urban area develops in a regular way following the master plan, this development can be regarded as sustainable because the social and economic driving forces play a major role in urban growth under the guidance of the master plan. By contrast, if urban growth is mainly manifested in a form of disorderly, outward sprawl and internal filling, it reflects that the master plan is out of operation and its socio-economic driving force is weak and loses support from government, organizations, and enterprises. In developing countries, the distribution and sprawling of slums provide an important indicator for assessing socio-economic sustainability. Slums are urban communities of low-and middle-income dwellers, with low, shabby and crowded houses, poor road system, and inadequate access to drinking water and sanitation (UN-HABITAT 2006). Slums are directly related to urban poverty and unemployment. The distribution of slums gives an index to the practice of urban planning and poverty reduction (Kohli et al. 2012). Urban poverty is manifested most conspicuously in the proliferation and expansion of slums (Kohli et al. 2016). Slum eradication has become the main goal of the UN millennium goals and the main content of the UN sustainable development 2030 goals (United_Nations 2015). Therefore, we can assess the socio-economic sustainability in urban area based on the distribution and evolution of slums.
The assessment of environmental sustainability using remote sensing technology is mostly reflected in the macroscopic conditions of human settlements. Here, we use NDVI to measure ecological health in forest area. NDVI is a vegetation index combining the red and near-infrared reflectance of a multispectral remote sensing sensor and efficiently characterize green vegetative cover and biomass in terrestrial ecosystem (Myneni et al. 1995;Lein 2014). When the red reflectance in the vegetation canopy is small and the near-infrared reflectance is large observed by remote sensing, NDVI increased to one in maximum, reflecting high vegetation coverage, high biomass and vigorous growth status, such as dense woodland. On the contrary, when the vegetation canopy has a large red reflectance and a small near-infrared reflectance, NDVI decreases towards zero, it reflects the decrease of vegetation coverage and biomass, and the poor growth status of plants, like sparse grassland or bare land. Long-term NDVI data can be used to assess land degradation (Chen et al. 2018a, b). Studies have shown that there is a good linear relationship between NDVI and precipitation in the growing season of grass land (He et al. 2015). Therefore, NDVI can be regarded as one of the important indicators for assessing and monitoring regional sustainability (Lein 2014). Located in the Margalla Hills national park in northwestern Islamabad, a large area of tropical and subtropical forests serves as a critical ecological barrier for the city, playing an important role in climate regulation, water conservation, soil erosion control, and biodiversity conservation. Thus, we can assess the ecological sustainability of Islamabad with the scale and health of forest area by using temporal NDVI data.

Results
In this section, we report the accuracy assessment results of image classification. We present variations and transitions among land cover types. Finally, we discuss the growing pattern of urban area in Islamabad.

Accuracy assessment of land cover classification
In previous study, we conducted a comparison among object-based post-classification change detection, backdating change detection associated with change vector analysis, and OBBCD. We used object-based image analysis to classify the Landsat images of the four periods, and then used the post-classification change detection results to assess the accuracy, showing that the overall accuracy of the three-time intervals was between 71.1% and 74.8%, and the kappa coefficient was between 0.66 and 0.71 (Method 1 in Table 3). Then we compare the results from integrating change vector analysis and object-based change detection with the fuzzy c-means-based result adopted in OBBCD. The total accuracy of the former is between 79.3% and 80.7%, and the kappa coefficient is between 0.74 and 0.77 (Method 2 in Table 3), while the total accuracy of the OBBCD approach is between 80.3% and 82.9%, and the kappa coefficient is between 0.77 and 0.81 (Method 3 in Table 3). The results show that the fuzzy c-means clustering method is the best for change detection.
Based on the OBBCD approach, we obtained land cover maps in 1990,1998,2009, and 2018 (see Fig. 2). The accuracy assessment showed that the results were highly available. The producer's accuracy (PA) and user's accuracy (UA) of the six land cover types in LD2018 were all larger than 87%. The overall accuracy and Kappa coefficient reached 90.3% and 0.882 (see Table 4). The overall accuracy and Kappa coefficient in LD1990, LD1998, and LD2009 were above 86.4% and 0.834. The PA and UA of the six land cover types were all larger than 85%. Among them, the PA and UA of impervious surface was above 85.4% and 83.2%.

Variations and transition of land cover types
These classification results performed regular variations of each land cover type from 1990 to 2018. We used three indices, including patch area, Percentage LANDscape (PLAND), and Largest Patch Index (LPI) from Fragstats (McGarigal and Ene 2015) to assess land cover changes. Other indices in Fragstats were not used because of their not-obvious variations in the time series. Figure 3 shows regular changes of the patch areas in the last 3 decades. Cropland (orange) increased from 275.21 km 2 in 1990 to 314.53 km 2 in 1998 and decreased afterwards to 227.75 km 2 in 2018. Accordingly, its PLAND increased from 29.1% in 1990 to 33.3% in 1998 and then decreased to 27.6% in 2009 and further decreased to 24.1% in 2018. Clearly, cropland area presented an overall downtrend. Its LPI increased from 7.46 in 1990 to 9.31 in 1998 and decreased to 6.36 in 2009 and to 2.69 in 2018, indicating an increased landscape fragmentation of cropland.
Woodland (green) demonstrated an increase 9.29% from 178.66 km 2 in 1990 to 195.26 km 2 in 2018. Its PLAND increased from 18.9% in 1990 to 19.7% in 1998, and remained the same 20.6% in 2009 and 2018, while its LPI changed small between 11 and 12. It implied that the fragmentation changed small. Most of increase happened in Zone-III.
Shrub and grassland (pink) decreased approximately in a half from 330.36 km 2 in 1990 to 163.58 km 2 in 2018. Its LPI decreased from 5.71 in 1990 to 3.71 in 1998 to 3.44 in 2009 and to 0.79 in 2018. Thus, the landscape fragmentation of shrub and grassland tended to increase.
Impervious surface (brown) demonstrated an increase of more than three times from 82.66 km 2 in showed a small increase from 6.49 km 2 in 1990 to 7.27 km 2 in 1998 to 7.52 km 2 in 2009 and to 7.64 km 2 in 2018. It took small PLAND values between 0.7 and 0.8 in the 3 decades. Further, the LPI of both bare land and water bodies varied between 0.21 and 0.69 and indicated a high landscape fragmentation.
It can be seen that, in Islamabad, urban growth leads to land cover change and induces the transition of cropland, shrub and grassland, and woodland towards built-up and paved area (Sudhira et al. 2004). In this study, the transition matrices of land cover types showed that the impervious surface was mainly converted from cropland, shrub and grassland, and bare land. Table 5 shows that the increment of impervious surface in 1998 were 4.73 km 2 from cropland, 12.63 km 2 from shrub and grassland, and 3.61 km 2 from bare land in 1990. Table 6 shows that the increment of impervious surface in 2009 were 23.91 km 2 from cropland, 28.7 km 2 from shrub and grassland, and 13.39 km 2 from bare land, respectively, in 1998. Table 7 demonstrates that the increment of impervious surface in 2018 was from 36.35 km 2 from cropland, 64.66 km 2 from shrub and grassland, and 15.37 km 2 from bare land, respectively, in 2009. These tables show that the transition among cropland, grass land, and bare land could be reversible, while the transitions towards impervious surface are irreversible. These transitions can be represented as Fig. 4.

Urban growth from 1990 to 2018
According to the two criteria for determining urban areas as stated in ''Extracting the urban area from remote sensing images'' section, we derived the urban areas in 1990, 1998, 2009, and 2018. In this extraction, the spatial adjacent relationship between Islamabad City and central Rawalpindi City was considered as well. It means that the impervious surface patches within 100 m of the boundary of Rawalpindi City refer to new urban areas of Islamabad City.
In the data in 1990 about impervious surface, there were a total of 1233 spatial patches (polygons), among which 88 meet the first criterion and 5 meet the second criterion. Thus, the urban area accounted for a total of 93 patches and reached 58.854 km 2 , taking 66.44% of the total impervious surface area and 6.22% of the total study area. As shown in Fig. 5, among the 88 patches that meet the first criterion, the largest unit had an area of 35.565 km 2 , covering the sector F-6, -7, -8, G-6, -7, -8, -9, and -10, and part of I-9, -10 adjacent to Rawalpindi, and a small part of F-5, G-5, H-8, H-9, and I-8 in Zone-I. The second largest area was in Defense Housing Authority (DHA) (it will be explained in next section) in Zone-V at the intersection of the highway E2 and N5, with an area of 3.079 km 2 . The third largest area was at the sector I-13 with an area of 2.409 km 2 around the boundary of both Zone-I and -II at the two sides of the highway N5. Five spatial patches which meet the second criterion were all adjacent to Rawalpindi City with areas less than 0.04 km 2 . The rest of the patches scattered in the five zones with areas between 0.1 km 2 and 1.0 km 2 .  From 1990 to 1998, the parts of urban growth contained 216 space patches with a total area of 23.576 km 2 . This accounts for 74.38% of the total impervious surface area and 8.71% of the total study area. This growth was mainly at the two sides of the highway N5 in both Zone-I and -II. The two largest patches were in the sector F-11, G-11 with an area of 1.354 km 2 and 0.978 km 2 , respectively. Additionally, the Bhara Kahu Town experienced a significant growth at the junction part of Zone-III and -IV.
From 1998 to 2009, urban growth involved 856 spatial patches with a total area of 81.924 km 2 . Urban area accounted for 88.05% of the total impervious surface area and 17.37% of Islamabad. They included (1) the sector  (2) Ali Bakhish Town near the border between Zone-I and -II; (3) Bhara Kahu near the border of Zone-III and -IV, and Ghori Town and Gulberge Green Town in Zone-IV adjacent to Rawalpindi City; and (4) DHA, MTH, Bahria Town along the highway N5 and E2 in Zone-V adjacent to Rawalpindi City. Between 2009 and 2018, urban growth contained 1326 patches with an area of 145.343 km 2 . It accounted for 93.71% of the total impervious surface area and 32.74% of the study area. The urban area in Zone-I and -II continued to expand, especially in the sector G-12, -13, -15, -16, H-9, -10, -11, -12, -13, and I-11, -12, -13 and -14. In Zone-II towns sprawled along the highway N5 at Tarnol Town. In Zone-IV, fringe sprawl along the eastern margin of Rawalpindi City made new settlements and slums between Ghori Town and Farash Town. In Zone-V, the urban expanded rapidly between N5 and E2. Prompt growth also happened in MTH, Bahria Town and DHA.
Our results show that the two criteria are effective to cover slums into the urban scope. With the increase of population and urban growth, slum area increased correspondingly in Islamabad. Before 1990, slums formed in the northeast Zone-I adjacent border of Zone-III. Among 58.85 km 2 of urban area, the slum area was 5.53 km 2 , accounting for 9.40% of the urban area. In 1998, the urban area was 80.47 km 2 . Of which the slum area reached 10.97 km 2 , increasing to 13.31% of the urban area. In 2009, the urban area reached 22.01 km 2 , while the slum area increased to 13.39% of the urban area, i.e., 22.01 km 2 . In 2018, the urban area reached 300.77 km 2 and the slum area increased to 28.98 km 2 and accounted for 9.36% of the urban area. Based on the statistics of urban area in the three periods, the growth rate of urban area increased at 5.01%, 9.04% and 9.83% per year. Correspondingly, the growth rate of slum area was slow down at 12.18%, 9.15%, and 3.52% per year. This implies that the urban is growing faster than slums. These results are showed in Fig. 5.

Discussion
At present, derivation of quantitative indices for sustainability assessment mostly rely on ground survey, census and socio-economic statistics, for which remote sensing technology only provides basic data, for example, urban spatial area. Further indicators from remote sensing oriented to sustainability assessment are still less reported. In this section, we discuss urban growth patterns and factors they affect. We assess the socio-economic and environmental sustainability of urban growth in Islamabad.

Importance of machine learning in object-based backdating change detection
The application of object-based backtracking change detection requires machine learning algorithms in many steps. Firstly, we need machine learning algorithm to classify segment objects from recent image.
In this study, random forest algorithm was proved to be more efficient, owing to sufficient sample data collection and abundant spectral and texture feature parameters of segment objects. Secondly, it is necessary to determine a classification threshold between change and non-change objects in the detection of object change. By comparison, we derived better results by means of fuzzy c-means clustering. Then we complete the type division of the change objects in the previous image using support vector machine. We found that this algorithm can achieve better results with limited training samples. Thus, machine learning plays a very important role in object-based backtracking change detection, and this role needs to be evaluated through the extensive use of various machine learning algorithms in the subsequent work.

Urban growth patterns
Urban growth usually includes a variety of patterns, such as fringe expansion, internal filling, and leapfrog development (Bhatta 2010). Different patterns imply different sustainability in socio-economic development. Following Bhatta (2010), the urban growth patterns in Islamabad can be divided into (1) village sprawl and leapfrog, (2) planning expansion, (3) fringe sprawl and infilling, and (4) merge.
The first type of growth patterns refers to villages sprawling outwards, upgrading to a town. Sprawl means a haphazard outgrowth of settlement areas with low population density, due to unrestricted, clumsy, and unplanned settlement growth . It is an inevitable pattern of settlement development under the background of active social and economic development. When a village far from the town grows and upgrades into a town through sprawling growth, it shows a leapfrog pattern of urban growth. This pattern appeared widespread around the sectors in Zone-I.
The second pattern of urban growth refers to urban expansion under the guidance of master plan. Made by the municipal government, master plan determines suitable use of urban land and functional zoning in urban space layout for a certain period. The plan affects the form and pattern in the planned area whenever the original landscape pattern is any kind of settlements or even non-urban land cover types, such as cropland, shrub and grassland, desert, or woodland. To increase this pattern, a road and highway network is built first. Then, the buildings and other complexes are constructed sequentially within the functional zoning areas, mainly happening in in Zone-I and -II.
The third pattern of urban growth involves sprawl from the urban fringe associated with infilling inside of a city. Fringe sprawl involves conversion of various land-use types at the urban fringe to buildings and roads. Infilling is inner sprawl associated with conversion from non-urban land-use types to urban and leads to increased building density. Fringe sprawl and infilling normally occur at the same time and present a macroscopic, disorderly pattern of buildings. It is usually accompanied by the emergence and spread of slums because of the lack of basic facilities, such as purified water supply, sanitation condition, and primary health center. This pattern presented in Zone-I, -II, and -IV around the downtown of Islamabad and north fringe of Rawalpindi City.
The fourth pattern is a combination of adjacent cities or towns, i.e. when they are close to each other and finally become one unit. The occurrence of this pattern usually involves merging two neighboring cities or towns of similar size into one, or expanding cities to the neighboring towns and villages in the periphery. The typical examples appeared in Zone-V.

Factors affecting urban growth
Behind these patterns, we recognize that there are mainly five factors pushing urban growth forward in Islamabad: population, transportation systems, master planning, industrial and real estate development, and neighbor urban effect.
The population in Islamabad has been growing rapidly in the past 30 years. According to the census, the population was less than a half of million until 1990, but it has reached two million in 2017 (Fig. 6). On this basis, the urban area has been steadily increasing. Thus, we can say that the increasing population pushes urban growth forward in Islamabad.
The second factor catalyzing urban growth is the transportation system, composed of roads and railways. Comparatively, there are relatively developed roads and railway systems in Zone-I, -II, -IV and the south part of Zone-V. They bring about more development opportunities than that in Zone-III. In Zone-IV and -V, towns expanded like chains along the Highway E2 and N5. Starting with Ghori Town, this chain included Ali Pur, Chatha Bukhtawar, Farash Town, and Bhara Kahu. One of railway lines leading to sector H-9 and I-9 helps people get in and out of the city and into jobs. The main railway line passes through Zone-V where Rawat Industrial Estate is located, carrying a huge quantity of raw material and cargo for manufactures.
The third factor promoting urban growth is the master plan. Zone-I and -II were developed following the master plan. The construction programs in the planning area provided new job opportunities and effectively supported urban growth. Many government agencies, parliaments, courts, and native and foreign organizations settled in Zone-I according to the master plan. To achieve the master plan, a road system was built to divide sectors with a regular grid system in the two zones, which supplies great convenience for improving the living and working conditions and attracted people settling down. Notably, Zone-III is not much developed because it was entrusted environmental protection in the name of Margalla Hills national park following the master plan.
The fourth factor is industrial and real estate development under the support of government and private enterprises. In 1988, the Pakistan Government announced industrial incentives to boost the economic growth in the country. One of the incentives was conducted at Potohar region in Zone-V. A joint venture initiated by both Export Processing Zone Authority (EPZA) and Rawalpindi Chamber of Commerce and Industries (RCCI) was established in 1980 and developed an export processing zone adjacent to the Rawat Industrial Estate in Zone-V. This is an autonomous organization with the mandate to plan,  1991, 1999, 2015violet columns to the urban area in 1990violet columns to the urban area in , 1998violet columns to the urban area in , 2009violet columns to the urban area in , and 2018 develop, manage, and operate export processing zones authorized by the ministry of industries. In the field of real estate, two private housing societies, the Defense Housing Authority (DHA) and Model Town Humak (MTH), have developed a large number of real estate projects and contribute many to urban growth in Zone-V.
The fifth factor is the neighboring urban effect. Rawalpindi is the contiguous city of Islamabad from the southeast and southwest side. Urban sprawling happened in Zone-I, -II, and Zone-IV and -V along the Highway E2 and N5 from the border of Rawalpindi. Several towns in Zone-IV sprawled from eastern fringe of Rawalpindi City. They passed through the Highway E2 towards eastern part of Islamabad, involving Burma Town, Ali Pur, Chatha Bukhtawar, Farash Town, Bani Galla, Bhara Kahu, Shehzad Town, and University Town.

Socio-economic sustainability
Urban construction in Islamabad is far from its target according to the master plan. Mostly because of its weak economic power, urban growth is uneven. The master plan aimed only at Zone-I and -II in the west part of Islamabad and omitted development planning in the rest zones. Even in Zone-I and -II, a lot of land remains to develop. Construction programs in no more than 20 sectors can be regarded as completed or close to completion. Further, urban economic development in Islamabad has been underpowered and urban development has been challenged by slum sprawl. Poverty and unemployment have long plagued the sustainable development of cities. Its poor population has also grown proportionately with the size of its slums.
The urban development was accompanied by slum sprawl in Islamabad. As early as the 1960s, indeed, when the capital city was founded in Islamabad, many slums appeared in the heart of Zone 1, such as the French Colony, Christian Colony, Musharraf Colony and settlements along streams and natural drains in sectors F-6, F-7, G-7, and G-8. Over the last 3 decades, all other zones ceaselessly arouse new slums except Zone-V, because the powerful development of industry and real estate in Zone-V do not reserve space for slums. In 1990, slums emerged in sectors E-11, F-6, F-7, F-11, G-7, G-8, H-13, I-13 and Nurpur Shahan in Zone-I; Tarnol, Katan, and Jhangisyedyan in Zone-II; Siadpur, Shah Allah Ditta in Zone-III; and Sohan, Tarlai Kalan, Ali Pur, Frash, Malpur, Lakhwal, Kirpa, and Shahpur along the highway E2 in Zone-IV (close to the eastern border of Rawalpindi city). This was because of easy to reach workplaces and saved transportation charges.
By 1998, new slums appeared in sector G-8, G-12, G-14, I-14, and Dhok pracha of Zone-I; Chatha Bukhtawar and Jhangi Sayedan in Zone-II; Bhara Kahu along E75 at the junction of Zone-III and -IV; Burma Town, Bhara Kahu, Chatta Bukhtawar, Kuri, Kothathyal, Jillani Town, Nogazi, and Partal Town in Zone-IV. From 1990 to 1998, many small slums also appeared at the two sides of the highway E75 with an area between 0.002 and 0.833 km 2 . The area developed into a large slum. Another prominent development was Model Town Humak (MTH) in Zone-V between the highway N5 and E2.
By 2009, new slums aroused in the sector E-13, F-12, G-12 in Zone-I. Original slums were filled in. Slums at Bhara Kahu, Burma Town, Sohan in Zone-IV merged from small scattered slums into a large concentrated slum. This period produced large slums extending along the E2, such as Khana Pull, Chatha Bukhtwar, and Bhara Kahu. In addition to the chronic lack of economic growth, another part of the flood of refugees caused by the war in the surrounding region. Islamabad is about 230 km to the border of Pakistan and Afghanistan. Since 2002, the United States has launched an anti-terrorism war in Afghanistan, carrying out a large-scale clearance of Al-Quaida (Taliban) forces. Series of border operations led to a large number of dwellings along border areas fled their homes and rushed into Islamabad and other cities. As a result, the area of slums in Islamabad doubled from 1998 to 2009 compared with the previous decade.
In 2018, new slums emerged in Zone-IV, including Pona Faqiran, Mochi Mohran, Morian. The former slums were infilled, sprawled merged, and thus contributed to the current pattern of urban growth. The problem of slum growth caused by the antiterrorism war remains.

Environmental sustainability of urban development
As seen from Fig. 7, the woodland shows a slow growth trend in the past 30 years. The main forest distribution in the north and East Mountains, a small number is located in lowland. Most of the newly grown woodland also extends outward or fills inward along the boundary of woodland areas. This can be partly attributed to the achievements made in regional ecological protection and construction, especially the active role played by the CDA and Margalla Hills Society.
Further, we used NDVI to check the health status in woodland. From the four Landsat images, we selected ten sites with a size of about 500 m 9 500 m in the extent of Margalla Hills national park. NDVI was calculated with red and near-infrared spectral reflectance from IM1990, IM1998, IM2009, and IM2018. It showed that NDVI tended to decrease in all these sites (Fig. 8). On average, it decreased from 0.668 to 0.551, reflecting a decrease trend of vegetation cover. It can be seen that although the forest area increased a small amount but the overall health became worse in the last 3 decades; so, its ecological quality tended to deteriorate. A research report shows that the temperature has shown a rising trend, while the precipitation has a large inter-annual change and is generally stable since 1980 (Khan and Fee 2014). The combination of the two factors leads to a regional trend of dryness, which may be the critical reason of the degradation of forest ecosystem as show above. The report also shows that underground water levels decline and tend to dry up, as urban demand for water increases. So environmental sustainability is under threaten because of climate change.

Conclusions
In this paper, we use the OBBCD method to obtain qualified land cover data of Islamabad from Landsat images acquired in 1990, 1998, 2009 and 2018. Research showed that over the last 3 decades, impervious surface area presented 273.10% increase from 88.579 km 2 in 1990 to 330.487 km 2 in 2018. Meanwhile, its impervious surfaces in Islamabad City raised from 9.36% to 34.94%. The growth is mainly from cropland, shrub and grassland, and bare land. We established two criteria for extracting urban land from impervious surface. Our study shows that the urban area increased 426.21% from 58.854 km 2 in 1990 to 309.697 km 2 in 2018. While the percentage of urban area in the study area increased from 6.22% to 32.74%. We recognized that machine learning plays a critical role in OBBCD.
On the basis of qualified urban land data, we evaluate the socio-economic sustainability and environmental sustainability of Islamabad City. We identified four urban growth patterns: village sprawl, planned expansion, fringe sprawl and infilling, and merging. We found that there are five factors pushing urban growth forward in Islamabad: population, transportation systems, master plan, industrial and real estate development, and neighbor urban effect. We recognized that the urban development in Islamabad is uneven. The urban growth in Zone-I and -II were mainly pushed by the master plan. The growth in Zone-IV is due to neighbor urban effect. The growth in Zone-V are based on industrial development and real estate. Although there was once a master plan in Zone-I and -II, the capacity of the city was insufficient, which leads to the development of only some sectors. Correspondingly, increasing slums spread in Islamabad with the rapid increase of population. The slum area increased from 5.534 km 2 in 1990 to 28.979 km 2 in 2018. The war in Afghanistan initiated by US troops led to a huge expansion of slums. Islamabad's socio-economic sustainability is perplexed by increasing slums.
The area and quality evaluation showed that woodland expanded 9.29% over the last 30 years because of ecological restoration programs in the Margalla Hills national park. Meanwhile, the quality of forest showed a decline trend. The mean NDVI decreased from 0.668 in 1990 to 0.551 in 2018, resulting in a deteriorated ecological quality.
Ethical approval All procedures performed in studies involving human participants were in accordance with the ethical standards of the institutional and/or national research committee. This article does not contain any studies with animals performed by any of the authors.
Informed consent Informed consent was obtained from all individual participants included in the study.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.