Enhanced normalized difference index for impervious surface area estimation at the plateau basin scale

Abstract. Impervious surface (IS) area is an important indicator of ecological environment condition in the basin. We propose an index for IS extraction [i.e., enhanced normalized difference impervious surfaces index (ENDISI)] by integrating the spectrum character of Landsat-operational land imager (OLI) images, and an automatic threshold selection method using the generalized Gaussian model. Dianchi and Erhai Basin are employed as study areas to test the ENDISI method at the plateau basin scale. The results show that: (1) the ENDISI can reduce the impacts of arid land, bare rock, and bare soil on IS extraction effectively; (2) ENDISI had a much higher separability degree between ISs and pervious surfaces compared with normalized difference built-up index, modified normalized difference IS index, and combinational biophysical composition index; and (3) the overall accuracy and kappa coefficient values of IS extraction via automatic threshold selection exceed 93.9% and 82.4%, respectively. Therefore, the ENDISI can serve as an effective index algorithm for rapid and high-precision IS extraction at the plateau basin scale.


Introduction
Impervious surfaces (ISs) have constantly expanded with significant effects on urban natural and human environments through rapid urbanization. IS cover is an important index used to quantify the influence of human activities and alterations of the environment. [1][2][3][4][5][6] The precise measurement of impervious surface area (ISA) helps control and improve the environment while also optimizing and regulating the scale of cities.
With the rapid development of sensor technologies and data processing techniques, remote sensing analyses and technologies have been continuously improved as more scholars have generated research results from IS extraction at different temporal and spatial scales. At present, many methods are used to extract IS cover such as spectral mixture analysis, [7][8][9][10][11][12] index methods [the normalized difference built-up index (NDBI), 13 the normalized difference impervious surface index (NDISI), 14 the modified normalized difference impervious surface index (MNDISI), 15,16 the normalized difference impervious index, 17 the combinational build-up *Address all correspondence to Kun Yang, E-mail: kmdcynu@163.com; Suozhong Chen, E-mail: 09153@njnu.edu.cn index (CBI), 18 the biophysical composition index (BCI), 19 combinational biophysical composition index (CBCI) 20 etc.], decision tree models, regression modeling, and other classification methods (maximum likelihood classification, object-oriented classification, artificial neural network classification, support vector machines, etc.). 21 The V-I-S model is a tool used for spectral mixture analysis, [22][23][24][25][26][27][28] which describes urban landscapes as a linear combination of three main components: ISs, vegetation, and soil. The V-I-S model has been widely used in many studies. 29,30 However, the method is affected by the selection of endmembers, rendering the results generated more uncertain. 31 Some index methods (e.g., the NDBI, CBI, and BCI) involve masking out water bodies or bare soil before extracting ISA. The modified normalized difference water index (MNDWI) is used to identify water bodies and to apply a water mask to remove water bodies from experimental results. Nevertheless, due to the diversity of bare soil components (various soil types, soil colors, and soil reflectance spectra), it is difficult to effectively identify soil and to set up a soil mask that removes it from the results. The NDISI uses the thermal infrared (TIR) band to account for the varying emissivity levels of ISs, which has clear advantages for urban mapping. However, as the TIR band has a low spatial resolution, extraction accuracy is limited by the presence of mixed pixels. Sun et al. 16 proposed an MNDISI that integrates the sharpened TIR band into the NDISI. When using the spectral discrimination index (SDI) as a performance evaluation index, the SDI values of MNDISI-derived maps are higher than those of the NDISI. Nevertheless, Chen et al. [32][33][34] found that the average temperature of ISs and bare soil is high across all seasons and that these surfaces always serve as urban heat sources. In addition, bright ISs are characterized by higher levels of surface reflectance than bare soil in the NIR and SWIR1 bands and in MNDWI values, and thus, the MNDISI-derived results are based on a portion of bare soil. Yang et al. 35 analyzed the spectral radiance of ISs and bare soil, obtained a normalized difference bare soil index (ρ) for the Dianchi basin, and facilitated the extraction of IS cover by eliminating most bare soils and vegetation using the INDBI algorithm. However, the choice of threshold has a stronger influence on IS extraction. The INDBI algorithm requires three thresholds to distinguish ISs from other underlying surfaces, rendering the results more uncertain. To highlight the four-major urban biophysical compositions (i.e., ISs, bare soil, vegetation, and water), a CBCI was proposed by Zhang. 19 Compared with NDBI, CBI, CBI, and NDVI, the CBCI has the significant correlation with IS and vegetation at different spatial resolution images, especially for high-resolution imagery. However, the CBCI cannot be applied to the region of Plateau Basin effectively, and there are still a lot of bare soil and arid land in the extraction results.
Previous studies have shown that bare soil can be easily confused with ISs. 36 Due to climate patterns, temperatures, humidity levels, rainfall levels, etc., there are not only many bare soil surfaces in plateau basins, but arid land (e.g., bare rock, gravel, and rocky desert) is also prevalent. Like bare soil surfaces, such surfaces are characterized by high surface reflectance levels and by similar spectral characteristics (especially in terms of NIR and SWIR1 values) to those of ISs, which seriously interferes with IS extraction.
This study addresses the impact of complex surface conditions on IS extraction for plateau basins. Based on previous research results, we analyze the spectral characteristics of an image and propose an enhanced normalized difference impervious surfaces index (ENDISI). The ENDISI can effectively eliminate the effects of arid land, bare rock, and bare soil on IS extraction. Using an automatic threshold selection method based on a generalized Gaussian (GG) model, the ISA are extracted from Landsat-8 imagery. Compared with NDBI, MNDISI, and CBCI, the results demonstrate that ENDISI has higher overall accuracy and SDI than other indices. Through case studies of the Dianchi and Erhai Basin in Yunnan Province, China, we estimate and analyze the IS extraction performance of the ENDISI at the plateau basin scale.

Study Area
The Dianchi and Erhai Basins were selected as study areas. The Dianchi Basin is positioned at 24°28′13″ to 25°27′26″ N and 102°29′14″ to 103°00′51″ E, covering a total area of 2843.16 km 2 and situated at elevations of 1652 to 2812 m [ Fig. 1(d)]. It is located within the Yunnan-Guizhou Plateau, which is mainly characterized by Karst plateau landscapes and red mountain original geomorphology. Annual average temperature and rainfall levels reach roughly 16.5°C and 1031 mm, respectively. The Dianchi Basin covers the Wuhua, Guandu, Xishan, Chenggong, and Panlong districts of the city of Kunming, which form the city's main areas of economic activity. Located southwest of Kunming, Dianchi Lake is the sixth-largest freshwater lake in China and the largest highland lake in Yunnan Province. The lake is an icon of the city of Kunming.
The Erhai Basin is positioned at 25°25′49″ to 26°25′55″ N and 99°49′56″ to 100°26′50″ E, covers a total area of 2608.73 km 2 , and is positioned at an elevation of 1847 to 4081 m [ Fig. 1(a)]. The climate of the Erhai Basin is that of a typical plateau continental climate with distinct wet and dry seasons, with a rainy season running from May to October and a

Data Source
As our main data source, we used Landsat operational land imager (OLI) images of the study area downloaded from the official U.S. geological survey (USGS) website. In the interest of obtaining quality analysis results, remote sensing images used were cloud-free or included few clouds. Their acquisition dates, spectral bands, spatial resolutions, and radiation calibration parameters are shown in Table 1.

Land use classification
Land use cover refers to the combination of natural foundations and surface elements covered by man-made structures, including vegetation, soil, lakes, wetlands, and various buildings with specific temporal and spatial attributes. In this study, we examine the geomorphic patterns of our study area and land uses are divided into eight categories: bright ISs, dark ISs, bare soil, arid land, forest, water, grass, and shaded surfaces. Here, the arid land class includes gravel land, bare rock, and stone desert land. The grass class includes urban green belts, lawns, and farmland with vegetation. Support vector machines (SVMs), known as maximum-margin classifiers based on the structural risk minimization principle of statistical learning theory, 38,39 are nonparametric statistical learning methods that have recently been used in numerous applications of image processing. 40 Mustafa Ustuner and B. DIXON's research 38,41 shows that SVMs can quickly accomplish high levels of classification accuracy with a small number of training datasets and are simple to apply. SVM was applied using environment for visualizing images (ENVI) software. Using ENVI version 5.1, we used a large number of training datasets for Landsat OLI images of a study area with eight categories of land use, and classification was performed using SVM classifiers.
To well obtain spectral information of different objects from images, we need to generate some spatial coverage samples 42 that are uniformly distributed throughout the study area. Webster and Oliver 43 provided a rule of thumb that 150 to 225 points might often be satisfactory for a variogram computed of isotropic variables in a normal distribution. In the SVMs classification results, the surface was divided into eight categories, so it needed at least 1200 to 1800 points. In addition, considering that some categories had fewer points and the cost saving. Therefore, 2000 samples were generated in this study, and they uniformly distributed in the study areas.
Taking the Dianchi Basin as an example. Using the "sp" and "spcosa" packages 42,44 under R version 3.4.4, the study area was divided into 2000 strata through multiple iterations (partitioning a spatial object into compact strata by means of k-mean), and the mean squared shortest distance of the strata was measured as 247868.2 m 2 . Corresponding to strata, we generated 2000 sampling points from the center and used attributes of the SVMs classification image, and thus, the 2000 sampling points were set as the land use prediction points for the study area. After checking the original and high-resolution images, the overall classification accuracy level was measured at 89.90% with 202 of 2000 prediction points being error points classified under incorrect categories (Table 2). We corrected the classification used for error points and obtained 2000 classification points. These can be considered verification points (Fig. 2).

Spectral curves of radiance and reflectance
Due to varying environmental conditions, the characteristics of radiated and reflected electromagnetic waves differ in wavelength. These are usually expressed as curves in a two-dimensional (2-D) geometric space. An abscissa denotes the wavelength (or band number), whereas an ordinate denotes the reflectance (or pixel value), which is called the spectral curve. Using the land-use classification data shown in Sec. 3.1.1, spectral curves of the radiance and surface reflectance were obtained from the images. Figure 3 shows the spectral curves for the Dianchi Basin.
Along the spectral curve shown in Fig. 3, the bright IS presents high values in each band, and waterbodies and shaded surfaces present relatively low values in each band. The spectral curve of the dark IS is similar to that of the bright IS, although the IS presents a correspondingly lower value.

Enhanced Normalized Difference Impervious Surfaces index
From the pixel spectral curves analysis shown in Sec. 3.1.2, we found differences between ISs and natural surfaces in partial band variations. First, bright and dark ISs present stronger reflectivity levels in the blue and coastal bands, whereas the reflectance values for water, vegetation, bare soil, arid land, and shaded surfaces are clearly lower than those of IS. Second, reflectance values for vegetation, bare soil, arid land, and shaded surfaces measured from the SWIR-1 band to the SWIR-2 band are quite low, whereas ISs maintain higher reflectance values with less change. Third, due to the strong effects of water vapor on the coastal band, we selected the blue band as an IS enhancement factor. To amplify the difference between ISs and natural surfaces, we selected the ratio of the SWIR-1 band to the SWIR-2 band and the MNDWI as the inhibitors. Finally, to obtain an index value of −1 to 1, the inhibited factor was multiplied by α. Thus, the ENDISI is written from Eq. (1): E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 3 6 7 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 3 0 8 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 2 6 1 where ρ Blue , ρ Green , ρ SWIR1 , and ρ SWIR2 are the surface reflectance values of the blue, green, first shortwave-infrared, and second shortwave-infrared bands, respectively. The ð•Þ Mean is the mean value of the image.

Automatic Threshold Selection Based on the Generalized Gaussian Model
The Kittler Illingworth (KI) threshold selection criterion assumes that observations are derived from a mixture of two normal distributions with respective means, variances, and proportions. 45,46 The criterion has been adapted to use GG distributions. 46 The optimization threshold T can be calculated from the cost function JðTÞ in the KI algorithm: JðTÞ ¼ þ HðΩ; TÞ − ½P PS ðTÞ ln α PS ðTÞ þ P IS ðTÞ ln α IS ðTÞ; E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 6 4 0 where hðX i Þ is the histogram of the log-ratio ENDISI with X i varying from L min to L max and with minimum and maximum values of the ENDISI incrementing at step l (l ¼ 0.01). Here, a j and b j are the constants related to shape parameters β j of the GG distribution for IS pixels and pervious surface (PS) pixels, respectively, and Γð•Þ is the Gamma function. The shape parameter β j (β j ≥ 0) determines the decay rate of the density function, and we used the corresponding estimation technique described in Refs. 47 and 48. Note that σ PS , σ IS , M PS , and M IS are the ENDISI standard deviations and the mean ENDISI of PS and IS pixels separated at T. P PS and P IS are the prior probabilities associated with PSs and ISs, respectively. HðΩ; TÞ denotes the entropy associated with the binary set of classes Ω ¼ fω PS ; ω IS g. HðΩ; TÞ is calculated as HðΩ; TÞ ¼ − where TÃ is determined by minimizing the cost function included in Eq. (4) with pixels presenting ENDISI values of greater than or equal to threshold (ENDISI ≥ TÃ) denoted as ISs and with those present values of less than the threshold (ENDISI < TÃ) denoted as PSs. In turn, an ISA map is extracted.

Other Impervious Surface Indices
To assess the performance of the ENDISI, three other indices, the NDBI, MNDISI, and CBCI were employed to conduct a comparative analysis. The equations of these three indices are as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 3 0 1 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 2 6 0 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 2 2 0 where ρ NIR and ρ SWIR1 are the surface reflectance values of the near-infrared and first shortwaveinfrared bands, respectively. Studies by Yang et al. 35 20 and the optimized soil adjusted vegetation index (OSAVI) was developed by Rondeaux et al. 49 The MBSI and OSAVI are calculated as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 1 6 ; 7 0 0 where ρ Red , ρ Green , and ρ NIR represent the surface reflectance values of the red, green, and NIR bands, respectively.

ENDISI Results
Using the ENDISI algorithm shown in Sec. 3 The ENDISI algorithm can suppress the effects of arid land, bare rock, and bare soil on IS extraction and can reveal information about ISs. The ENDISI results show that the values for waterbodies, shaded surfaces, and forest are very low, most of which fall below −0.25. Bright and dark ISA present higher ENDISI values in plateau basins, and they can be separated from bare soil and arid land areas. It is easy to observe the spatial distribution of ENDISI, with high ENDISI values concentrated in major urban areas and at airports, whereas the low ENDISI values mainly found throughout plateau lakes (Dianchi Lake and Erhai Lake) and along the borders of basins. We must add that as the reflectance of snow is very high, the numerous snowy surfaces

Statistical Analysis
To further examine the ENDISI's capacity to distinguish between land covers, we calculated the ENDISI values for the 2000 land use classification points corrected in Sec. 3.1.1. A statistical analysis of the ENDISI value distributions found for various land use categories was in turn conducted (Fig. 6).
The statistical results show that bright and dark ISs present high values, whereas waterbodies, shaded surfaces, forested land, and arid areas present lower values. Most bare soil areas present lower ENDISI values than ISs, and only a small portion of bare soil areas present ENDISI values that overlap with those of dark ISs. In addition, some farmland areas present ENDISI values that overlaps with those of IS areas. Local farmers have built many translucent artificial greenhouses on farmland, which has enhanced surface reflectance levels and thus affected ENDISI values.

Automatic Threshold Selection
Threshold selection is central to ISA mapping, as it directly affects the accuracy of ISA mapping. Previous methods involve manually selecting a threshold through continuous debugging. This is not only time consuming but is also affected by human interference, often causing the thresholds of different images to vary. To address the above problems, we adopt an automatic threshold selection method that uses the minimum error threshold method as a constraint to enable automatic, rapid, and efficient selection over a relatively optimal threshold. Figure 7 shows the automatic threshold selection process used for ENDISIs of the Dianchi [ Fig. 7(a)] and Erhai Basins [ Fig. 7(b)]. The KI threshold selection criterion based on the GG model proved to be effective in this experiment. Automatic thresholds applied in the ENDISI histogram of the Dianchi and Erhai Basins are measured at −0.16 and −0.09, respectively.

Separation of ISs and PSs
We reclassified the 2000 classification points corrected in Sec. 3.1.1. The classification points are divided into IS and PS categories and are used as verification samples. The ISs contain bright ISs and dark ISs, and the PSs include natural surfaces such as bare soil, arid land, forested land, waterbodies, grassland, and shaded surfaces (as shown in Table 3). Figure 8 shows ENDISI histograms for IS and PS samples for the Dianchi [ Fig. 8(a)] and Erhai Basins [ Fig. 8(b)]. Visually, ISs and PSs present respective distributions and GG distributions of two different shape parameters. Although some overlapping ENDISI values between IS and PS are observed, the ENDISI separates ISs and PSs quite well.
The spectral discrimination index (SDI) can distinguish between two groups of data and can quantify a histogram's separation with higher SDI values denoting more separable results. 50 For more detailed information on the SDI, please see Sun et al. 16 and Pereira et al. 51 From the verification samples, ENDISI SDI values were calculated via ISA extraction. ENDISI SDI values from the Dianchi and Erhai Basins were measured at 1.499 and 1.783, respectively.

Overall accuracy
From our verification samples, an accuracy assessment was conducted using the error matrix approach. 52 Here, the overall accuracy and overall kappa coefficient of the ENDISI were calculated for any given threshold (step ¼ 0.01). For the Dianchi Basin, the optimal threshold is −0.16 with 93.95% overall accuracy and 82.47% kappa coefficient. For the Erhai Basin, the optimal threshold is −0.07 with an overall accuracy level and kappa coefficient of 98.90% and 93.21%, respectively. At each given threshold, the overall accuracy and kappa coefficient values are shown in Fig. 9.

Comparative Analyses with Other Indices
Taking the Dianchi Basin as an example, NDBI (replacing ρ SWIR1 with ρ SWIR2 in the NDBI equation), MNDISI, and CBCI images were obtained using the equations shown in Sec. 3.4. Figure 10 shows the distribution of three indices for different land covers. Bare soil and arid land surfaces are the main interference sources of the three algorithms. Especially in MNDISI, the average values of bare soil and arid land surfaces are larger than bright IS and dark IS, which directly affect the mapping efficiency of MNDISI. Compared with MNDISI, NDBI, and CBCI can better resist the impact of vegetation information (forest and grass surfaces). The SDI value indicates the separability between ISs and pervious surface. Among the four ISA extraction algorithms, ENDISI has the highest overall accuracy and SDI value (Table 4), reaching 93.95% and 1.499, respectively. Moreover, the relatively high SDI value (0.906) of the NDBI was also lower than that of the ENDISI, which was supported by the small overlaps of the histograms between the ISs and pervious surfaces [ Fig. 11(a)]. Additionally, the histograms of  MNDISI between the ISs and pervious surfaces displayed obvious overlaps [ Fig. 11(b)], which was supported by a relatively low score of the SDI value (0.498). As shown in Fig. 10(c), some water bodies present CBCI values that overlap with those of IS areas, which was supported by a relatively low SDI value (0.705).  Fig. 13(c)], almost all pixels are classified as ISs with the exception of a few vegetation-covered areas. Additionally, there are some large areas of sand factories in Fig. 12(c), lots of sand and gravel have been accumulated on the ground. Because they are the main raw materials for ISs, the factory areas are also classified as ISs.

ENDISI for Landsa-TM/ETM Imagery
The ENDISI can be well applied to Landsat-TM or ETM images. With the automatic threshold extraction algorithm in Sec. 3.3, we can get the ISA distribution map quickly. Meanwhile, Figs. 14 and 15 were the consequences of ENDISI extraction with Landsat TM and ETM images as data sources, respectively.      The experimental results illustrate that the ENDISI can well distinguish the bare soil and the arid land from the IS in the Landsat TM and ETM images, which overcomes the effect of the bare soil and arid land elegantly [Figs. 14(a) and 15(a)]. For one thing, in high ISA, the ENDISI can well differentiate buildings, airports, streets, and other man-made structures from natural surfaces [Figs. 14(b) and 15(b)]. For another, in low IS density areas, the ENDISI can extract the main roads and villages from the images [Figs. 14(c) and 15(c)]. However, owing to low spatial resolution and mixed pixel images, some narrow roads and mud roads in the study area were divided into pervious surfaces [ Figs. 14(a) and 15(a)]. Generally speaking, the ENDISI can achieve good IS extraction effect on Landsat series images.

Discussion
With accelerated urbanization, ISA have rapidly expanded, representing a major source of environmental pollution in plateau basins. In this study, we extracted ISs of plateau basins from the ENDISI index and via automatic threshold selection. Through our experimental analysis (Sec. 4.5), we found that bare soil, bare rock, and arid land areas interfere with IS extraction, as their spectral curves are similar to those of ISs and especially to those of bands 5 (NIR) and 6 (SWIR-1). Therefore, ISA extracted via the normalized difference build-up index (NDBI) method include large areas of bare soil, bare rock, and arid land. In addition, the average surface temperature and brightness temperatures of bare soil are higher than those of ISs in the Dianchi Basin, 32 and the surface temperatures of the arid land are also similar to those of ISs. Therefore, the results of the MNDISI and NDISI algorithms also show large expanses of bare soil and arid land, not suitable for IS extraction in plateau basins. Furthermore, the reflectance spectral curves of different soil types for different regions are discrepant. Therefore, the improvement of MBSI can provide a better performance for CBCI.
The GG distribution is a symmetric distribution, which uses the δ function and the uniform distribution as a limit form. The GG distribution can fit statistical data well and has a wide range of applications once its shape parameters are regulated. 53,54 In this study, ISs and PSs were found to be well expressed by the GG distribution of a plateau basin, and it was found that an automatic threshold can be derived from the minimum error threshold method described in Sec. 3.3. In turn, automatic ISA mapping can be performed. This method works well when applied to Landsat imagery of a plateau basin, and the feasibility of using the ENDISI to extract ISs from other remote sensing images will be tested in the future studies.
IS extraction using the ENDISI algorithm remains limited in certain respects. During urban construction processes, attention is paid to the harmonious construction of green ecology. Almost all pixels of a major urban area are mixed pixels, creating interference with IS extraction. Urban green belts and grassland are often positioned adjacent to the dark ISs, and some agricultural areas adopt greenhouse structures for crop growth. The plastic roofs of greenhouses are composed of semi-transparent artificial materials, altering the reflectivity of vegetable farms. Therefore, some grassland areas in the Dianchi Basin present higher ENDISI values [as shown in Fig. 8(a) in which the histogram distributions of IS pixels are overlain with PS pixels]. Experimental results and analysis show that ENDISI algorithm can effectively reduce the impact of soil and arid lands. Therefore, the ENDISI can not only serve the ISA mapping of plateau basin but also be well applied to the ISs extraction of other urban areas.

Conclusions
This paper proposes an ENDISI that can efficiently extract ISs of a plateau basin. We adopt an automatic threshold selection method based on the GG model to quickly and accurately select a threshold for IS extraction from an ENDISI image.
The quantitative analyses show that the ENDISI performs relatively well in separating ISs from PSs, generating SDI values for the Dianchi and Erhai Basins of 1.499 and 1.783, respectively. Compared with NDBI, MNDISI, and CBCI, the results show that ENDISI has a much higher separability degree between ISs and pervious surfaces. The automatic threshold method can be used to quickly and automatically select a threshold that is similar to the optimal threshold with an accuracy deviation of <0.4%, thus ensuring IS extraction. Automatic threshold selection eliminates uncertainties associated with manually selecting thresholds for efficient and effective IS extraction. When using the automatic threshold, the overall accuracy and kappa coefficient levels derived from IS extraction from the study areas exceed 93.9% and 82.4%, respectively.
The ENDISI and automatic threshold selection method proposed in this paper can be used as effective tools for the efficient and automatic extraction of ISs of plateau basins and are suited to the large-scale monitoring of IS changes.