Using long time series of Landsat data to monitor impervious surface dynamics: a case study in the Zhoushan Islands

Abstract Islands are an important part of the marine ecosystem. Increasing impervious surfaces in the Zhoushan Islands due to new development and increased population have an ecological impact on the runoff and water quality. Based on time-series classification and the complement of vegetation fraction in urban regions, Landsat thematic mapper and other high-resolution satellite images were applied to monitor the dynamics of impervious surface area (ISA) in the Zhoushan Islands from 1986 to 2011. Landsat-derived ISA results were validated by the high-resolution Worldview-2 and aerial photographs. The validation shows that mean relative errors of these ISA maps are < 15   % . The results reveal that the ISA in the Zhoushan Islands increased from 19.2     km 2 in 1986 to 86.5     km 2 in 2011, and the period from 2006 to 2011 had the fastest expansion rate of 5.59     km 2 per year. The major land conversions to high densities of ISA were from the tidal zone and arable lands. The expansions of ISA were unevenly distributed and most of them were located along the periphery of these islands. Time-series maps revealed that ISA expansions happened continuously over the last 25 years. Our analysis indicated that the policy and the topography were the dominant factors controlling the spatial patterns of ISA and its expansions in the Zhoushan Islands. With continuous urbanization processes, the rapid ISA expansions may not be stopped in the near feature.


Introduction
Coastal islands are located in a transitional environment where land and ocean interact. Such an environment is vulnerable to human activities due to limited natural resources. Under urbanization and intensive exploitation activities, the coastal zone ecosystem health in Zhejiang province is facing serious threats. 1 During the past 20 years, large areas of natural landscape of the Zhoushan Islands have been increasingly replaced by impervious surfaces. Impervious surface area (ISA) can be defined as any material that prevents the infiltration of water into the soil, such as paved roads, driveways, parking lots, buildings, rooftops, and so on. 2 The expansion of impervious surfaces can affect the hydrological cycle, water quality, local climate, and biodiversity. 3 It was found that most stream quality declines if >10% area of a watershed is impervious. 4 Quantification of impervious surfaces in landscape has become increasingly important with growing concern of its impact on the environment. [5][6][7][8] Remote sensing has long been used for land use/land cover (LULC) classification, and it allows up-to-date, spatially explicit estimates of urban surface imperviousness over large areas. [9][10][11] Ridd proposed the classic vegetation-impervious surface-soil model to parameterize biophysical composition of urban environments. 12 Carlson and Arthur retrieved the fraction impervious surface of Chester Country based on the inverse relationship between impervious surfaces and vegetation. 2 Wu and Murray estimated impervious surface distribution of the metropolitan area of Columbus, Ohio, through a fully constrained linear spectral mixture model. 13 Xian and Crane determined impervious surfaces of the Tampa Bay watershed in Florida, to assess urban sprawling, using a regression tree algorithm. 14 Object-based image analysis has been increasingly used to map impervious surfaces due to the advent of high-resolution satellite imagery. [15][16][17] Xu proposed a novel index, normalized difference impervious surface index (NDISI), to estimate impervious surfaces of Fuzhou City with a high degree of automation compared to the majority of the present methods. 18 To sum up, the methods for mapping ISA with remote sensing technology can be classified into four groups: visual interpretation or multispectral classification, complement of vegetation fraction, spectral mixture analysis, and band combination method. Previously mentioned studies were mostly applied to inland regions and they might not be suitable for islands. For islands, in addition to types of woodland, arable land, urban or built-up areas, rural settlements, and water bodies, there are also brine pan, tidal zone, and aquaculture, the spectral curves of which can easily be confused with impervious surfaces under the influence of water. 19,20 Zhou and Wang developed a multiple agent segmentation and classification algorithm that included the steps of a multiple agent segmentation, shadow-effect, multivariate analysis of variance-based classification, and postclassification to extract ISA maps in the state of Rhode Island, USA, using Quickbird-2 imagery. Besides, many researchers have done some work in island land-cover analysis using remote sensing and geographic information system. [21][22][23][24][25] Chen et al. explored the temporal composition and spatial configuration of the land-cover trajectories in the ZS and its surrounding islands from 1986 to 2000, and revealed that the land cover had changed rapidly in the ZS. 23 However, little attention was paid to the evolution of island impervious surface over time, and long-term impervious surface mapping and change analysis in islands have not been conducted. Therefore, this study will take the Zhoushan Islands as a case study to map and analyze impervious surface dynamics over the past 25 years.
Choosing one or more end members to represent multiple impervious surfaces is not easy. This is particularly true for applying low-and medium-resolution remote sensing data. 26 The construction of the NDISI needs thermal data, and therefore the remotely sensed data without thermal bands are useless in the indicator. Visual interpretation or multispectral classification is mainly applied to high-resolution imagery. Because of the inverse correlation between ISA and vegetation cover in urban areas, one potential approach for impervious surface extraction is through vegetation distribution. 2,[27][28][29] The methods for acquiring vegetation fractions from remote sensed data are quite mature following many years large-scale experiments. [27][28][29][30][31] Here, we assume that within the urban (or developed) pixels, all areas without vegetation covers belong to ISA and ISA are only distributed in these urban (or developed) pixels. Therefore, mapping ISA needs two data sources: vegetation cover and urban distributions. Based on Chen's study, the newly developed areas in 1986,1990,1995,2000,2006, and 2011 can be acquired. 23 Another parameter, fractional vegetation distribution, can be calculated by the scaled normalized difference vegetation index (NDVI) or tasseled cap transformation. 28,32-34 Therefore, we decided to combine the time-series classification and fractional vegetation cover to map and monitor impervious surface dynamics of the Zhoushan Islands.
The objectives of this research are (1) to build time-series ISA maps in the study area, (2) to estimate ISA and its changes over different time periods, (3) to examine conversions between impervious surfaces and other LULC types, and (4) to analyze temporal and spatial variations of ISA distribution in the new urban areas.

Study Area
The study area is part of Dinghai District, Putuo District, and Daishan Country (29°55'-30°12'N, 121°48'-122°20'E) of Zhoushan archipelago new area, located in the northeast of Zhejiang province, China. It includes ZS and its surrounding inhabited islands, namely Jintang (JT), Xiushan (XS), Cezi (CZ), Changbai (CB), Changzhi (CHZ), Damao (DM), Aoshan, Xiaogan-mazhi (XG), Panzhi (PZ), and Lujiazhi (LJZ), covering an area of ∼703.34 km 2 (Fig. 1). The basic information of these islands is listed in Table 1. The Zhoushan Islands have a hilly landscape and their climate belongs to the monsoon-influenced subtropical marine system. The mean annual temperature is around 17°C, and the average annual rainfall is 1300 to 1500 mm. Coniferous forests, evergreen broad-leaved forests, and deciduous broad-leaved forests are widely distributed on this archipelago. The capital, Zhoushan city, has been one of the pioneering cities in the Yangze River Delta area since 1988. With the Zhoushan cross-sea bridge between Zhoushan and Ningbo city, the study area is subject to rapid changes in transportation, ports, tourism, and fisheries. Economic prosperity in the Zhoushan Islands has generated a significant competition over limited land resources. Quantifying the ISA and analyzing impervious surface change during the past 25 years in the Zhoushan Islands are therefore of great importance to characterize the influence of anthropogenic activities in this area.

Data and Preprocessing
Remotely sensed data sets used in this study are listed in Table 2. The selected images were georeferenced into a common transverse mercator projection with geometric errors <1 pixel. Then the digital numbers of Landsat images were retrieved to apparent reflectance by using the model of Chander et al. 35 The nearest-neighbor resampling technique was used to resample the images into a pixel size of 30 m by 30 m for Landsat images, 5 m by 5 m for SPOT5 multispectral images, and 0.5 m by 0.5 m for Worldview-2 multispectral images. We used Pansharp fusion method from ENVI EX 4.8 to merge images of SPOT5 and Worldview-2.
CORONA images from the year 1970 were used to map early historical land cover over the study region. Landsat TM∕ETMþ imagery and SPOT5 were used to quantify the LULC and    (Table 3). Other data sets used in this study include differential global positioning system data and the land-use vector data from Zhejiang Administration of Surveying Mapping and Geoinformation, China; these are utilized as the reference data for classification purposes.

Estimating ISA
ISA is generated by applying the following formulation. 2,36 As mentioned earlier, we assume that ISAs are only distributed in pixels being classified as "developed" and it is the fraction of each pixel that is not covered by vegetation.
where the subscript dev indicates that the quantity is defined only for pixels classified as "developed" (i.e., urban/suburban) and F r indicates fractional vegetation cover in a pixel.

Fractional vegetation cover
To obtain F r , we first compute NDVI. NDVI is defined as where ρ r is reflectance in the red band (630 to 690 nm) and ρ nir is reflectance in the near-infrared band (760 to 900 nm). To remove biases from atmospheric influences on NDVI, we used two references, one for bare soil (NDVI s ) and one for full vegetation (NDVI v ), respectively, to rescale the calculated NDVI. 33 The rescaled vegetation index (N Ã ) is calculated as In this research, we overlaid NDVI maps for each time on the corresponding classification map to derive NDVI v in the land type of "woodland" and NDVI s in the land type of "dryland," respectively. Time-series classification maps are discussed in Sec. 3.1.2.
Field experiments and observations indicated that F r has a linear relationship with N Ã2 (Refs. 32, 36, and 37) and we used this equation to quantify F r from N Ã .  23,25 The study area was classified into eight classes, including woodland, arable land (i.e., dry land and paddy field), built-up areas, rural settlements, water body, aquaculture, brine pan and tidal zone. These land-use types and interpretation criteria for Landsat images were used for classifying historical LULC from CORONA image in 1970. The 2005 classification map was directly from the vector data in 2005. The new LULC map for 1990 was generated with the same method as Chen et al. 23 The samples for supervised classification were selected from pixels without change by comparing the LULC maps before and after 1990, respectively. Using the land-use vector map of the year 2005 as a reference, Landsat TM and SPOT5 images were classified for land use in 2006 and 2011, respectively. The classification accuracies were higher than the standard of 85% stipulated by the United States Geological Survey classification scheme. 38 For lack of high-resolution data in 1986 and 1995, classification accuracies in these two periods were not quantified.
Here the urban/built-up areas and rural settlements were chosen to generate masks of developed pixels. The nondeveloped areas were reclassified as 0% impervious surface. Moreover, modified normalized difference water index (MNDWI) 39 was combined to exclude some rivers ignored in developed areas as follows. Practice has proved that the pixels are regarded as water when satisfying 0.18 > MNDWI > 0.1.
where ρ g and ρ m1 are the apparent reflectance of TM∕ETMþ green band (520 to 600 nm) and middle-infrared band (1550 to 1750 nm), respectively.

Accuracy Assessment
The accuracy assessment involved the following steps. (1) Samples were randomly chosen from the same period's high-resolution data, but excluding the shadow pixels and pixels that are totally different from low-resolution images.
(2) Classifications on the sampled pixels of high-resolution images were conducted with integrated methods including visual interpretation, maximum likelihood classification, and vector quantization. The overall accuracy and the kappa coefficient of these classifications on samples should be >90% and 0.89, respectively. (3) All 1-m pixels classified as impervious surfaces were tallied within TM pixels to compute ISA. (4) Accuracy assessment indicators of the mean relative error (MRE) and Pearson correlation coefficient (R) were selected. 40 Table 3). The accuracies on ISA mapping are listed in Table 3. On the whole, the correlation coefficient is >0.6 and MRE is <15%. There are some uncertainties on these evaluations on the ISA mapping. The different image acquisition time between reference high-resolution images and TM∕ETMþ images may show a spectral difference of the same object, which produces the different ISA, though in the same invariant locations. Besides, the assumption that the impervious surface pixels would not return to their pervious surface type at a later time may not be true spanning such a long period.

Differences in ISA Changes Among Major Islands
The percentages of ISA distribution for the entire study area from 1986 to 2011 are shown in Fig. 2. The total ISA in the Zhoushan Islands increased from 19.2 km 2 in 1986 to 86.5 km 2 in 2011. The growing rate of ISA is continuously increasing during 1986 to 2011 (Fig. 3) and the last period has the highest growing rate of 5.59 km 2 per year.
Because of the difference in location, geomorphology, and population among these islands, they have different changes in total ISA and the percentage of ISA. ZS has the largest increase in total ISA, followed by JT, XS, CHZ, CZ, XG, LJZ, PZ, and, DM (Table 4). In terms of change in percentage of ISA, LJZ goes to the first place, followed by CHZ, XG, CZ, XS, PZ, ZS, JT, and DM (Table 4). Islands closer to ZS normally have a larger increase in ISA than farther islands. For example, CHZ, which is 0.3 km from ZS, has seven times more increase in the percentage of ISA than DM, which has the same size as CHZ but is located far away from ZS. Some smaller islands, such as CHZ, XG, and LJZ, had higher increases in the percentage of ISA than the main island (ZS), even though the policy of "construct the major island and to emigrate from minor islands" was started in 2000. Because of the limited plain area in ZS, the land requirements for industries in this region will fall on nearby islands with ideal landscape and locations.

Land Conversions to ISA
The land-cover change analysis from 1970 to 2011 has been conducted in previous research. 25 Here we focused on the conversions from other land covers to ISA in the Zhoushan Islands. The LULC maps in the study area from 1986 to 2011 are shown in Fig. 4. The most notable change was built-up areas, which grew by twice the area in 1986, followed by arable land and tidal zone, which shrank by 40.27 and 30.05 km 2 , respectively. Aquaculture totally increased 10.52 km 2 , and brine pan decreased 8.57 km 2 from 1986. The major land conversions in the Zhoushan islands over the last 25 years were from tidal zone to aquaculture and brine pan, from brine pan to built-up areas and aquaculture, and from arable land to built-up areas. Overall, more and more natural landscapes were disturbed by human activities and the land use had been intensified.
The percentage of ISA was classified into three ranks: 0 to 60%, 61 to 80%, and 81 to 100%. By integrating ISA maps, high-resolution images, and land-use investigation maps, we found that high-percentage ISA (81 to 100%) mainly covered industrial land, port land, and some urban high-density buildings; middle-percentage ISA (61 to 80%) involved residential areas outside the city center, hydraulic construction sites, and roads; and low/lower percentage ISA (<60%) was composed of green parks and suburban woodlands. During the past 25 years, Zhoushan has experienced a rapid urbanization. The population of urban areas has doubled since 1986. 42 The increasing population in Dinghai District and Putuo District drove the conversions of natural and arable lands into built-up and urban areas, including infrastructures and facilities, e.g., hospitals, roads, and residential buildings. It generated some new middle-percentage ISA areas in the suburban areas and some new high-percentage ISA areas in the inner city. Besides urbanizations near the big city, there are a growing number of high-percentage ISA areas along the shoreline, e.g., ports, piers, and industrial lands. These were largely converted from tidal zones.   However, for the periods 1990 to 1995 and 1995 to 2000, large areas of new impervious surfaces formed within and around existing urban areas. After 2000, it was observed that the impervious surface of Dinghai District grew from the inner city to the west and southeast, and gradually merged with the Putuo District. From 2006 to 2011 the impervious surfaces in ZS, XS, CHZ, CZ, and JT islands experienced a rapid growth along the shore. The driving forces of these spatial and temporal patterns of ISA were related to policies and geography. First, the policy that "construct the major island and to emigrate from minor islands" caused the increase in impervious surfaces in ZS from 1990 to 2005. During this period, increasing populations required a growth of infrastructures, including energy, hospitals, roads, subsidized housing, as well as village-owned enterprises in towns. Second, the strategy of "to thrive the city by relying on the port" in 2005 boosted the exploitation of ports, piers, and industrial lands along the shoreline in this region. Finally, the topography also had a certain impact on the distribution of ISA. In ZS, there are hills and mountains in the middle and only small areas of plains distributed in the periphery [ Fig. 1(b)]. The plain land existing in the nearby islands gave an expansion chance for Zhoushan urban areas, and therefore, CHZ and LJZ are functioning as urban extensions of Dinghai District and Putuo District, respectively. With the expansions of urban area, the frequencies of the percentage of ISA (ISA%) were changing over time (Fig. 6). Normally, the mean and mode (or the peak) of the ISA% is becoming higher and higher since the year 1986 (Table 5).   (1986 to 1990, 1990 to 1995, 1995 to 2000, 2000 to 2006, 2006 to 2011) in the new urban areas. The ISA of 1986 was conducted as the reference. Fig. 6 The frequency distribution of ISA in the increased urban area, 1986 to 2011.

Conclusions and Future Works
This study generated ISA maps of the Zhoushan Islands using multiple satellite data from 1986 to 2011. These maps show that the total ISAs in the Zhoushan Islands expanded from 19.2 km 2 in 1986 to 86.5 km 2 in 2011, at an average increase rate of 2.70 km 2 per year. The increasing rate of ISA in these islands had large spatial heterogeneities because of differences in the topography, size, and location. The dominant land conversion in these islands was from tidal zones and arable lands to urban areas with middle and high percentages of ISA. The expansions of ISA occurred unevenly in spatial dimension and most of the new ISA took place in peripheries of these islands. The mean and mode of ISA in urban area increased gradually from 1986 to 2011. The policy and the topography may contribute to the pattern of ISA and its expansions for different periods. Results from this research can be used by policymakers for island management, planning, and for ecological and hydrological modeling to determine the effect of the increasing ISA on coastal environments of the Zhoushan Islands. Several issues are not covered in this study. One issue is the use of spectral mixture analysis classification approach. Although we have estimated the ISA and its change rate for the Zhoushan Islands from 1986 to 2011, and have captured the general impervious surface dynamic using middle-resolution satellite data, the use of different spectral endmembers for tidal zone, brine pan, aquaculture, and impervious surfaces may improve the classifications on ISA. Another issue is the effects of spatial resolution of remotely sensed data on estimating island ISA. The high-resolution data may better detect the imperviousness of the ring belts, from the central to the urban-rural peripheral, than the Landsat data. 26 In this study, the 2011 SPOT5/ Worldview-2 data were only used to validate the extraction of impervious surfaces. Highresolution satellites might be used in the near future when they become economically affordable for large-scale mappings.