A Normalized Urban Areas Composite Index (NUACI) Based on Combination of DMSP-OLS and MODIS for Mapping Impervious Surface Area

Abstract

Mapping Impervious Surface Area (ISA) at regional and global scales has attracted increasing interest. DMSP-OLS nighttime light (NTL) data have proven to be successful for mapping urban land in large areas. However, the well-documented issues of pixel blooming and saturation limit the ability of DMSP-OLS data to provide accurate urban information. In this paper, a multi-source composition index is proposed to overcome the limitations of extracting urban land using only the NTL data. We combined three data sources (i.e., DMSP-OLS, MODSI EVI and NDWI) to generate a new index called the Normalized Urban Areas Composite Index (NUACI). This index aims to quickly map impervious surface area at regional and global scales. Experimental results indicate that NUACI has the ability to reduce the pixel saturation of NTL and eliminate the blooming effect. With the reference data derived from Landsat TM/ETM+, regression models based on normalized DMSP-OLS, Human Settlement Index (HSI), vegetation adjusted NTL urban index (VANUI), and NUACI are then established to estimate ISA. Our assessments reveal that the NUACI-based regression model yields the highest performance. The NUACI-based regression models were then used to map ISA for China for the years 2000, 2005 and 2010 (Free download link for the ISA products can be found at the end of this paper).

1. Introduction

Over the past 30 years, global urban population has increased from 1.5 billion in 1975 to 3.4 billion in 2010. This growth trend is predicted to continue into future decades according to the United Nations. Global urban population will increase to 6.3 billion by 2050; nearly double the current urban population [1]. The massive immigration to cities has resulted in correspondingly rapid urban expansion. Although urban settlements occupy only a small fraction of the world’s surface, they have a significant effect on climate, biogeochemical cycles, biodiversity and hydrology at local, regional and even global scales [2]. Unprecedented urban expansion has brought about a series of environmental and ecological problems such as encroachment on agricultural land, destruction of ecosystems, water shortage, air pollution and urban heat islands [3,4,5,6]. Urban land is increasingly being recognized as an important yet poorly quantified component in global change models [7]. Accurate information about urban dynamics is crucial for revealing the relationships between urbanization and its environmental consequences [8,9].

Remote sensing techniques, which can provide quick, periodic revisits and large area coverage, have proven to be a useful tool for mapping urban land [10]. Satellite images have been extensively used in monitoring urban dynamics [11,12,13]. In most previous studies, urban area mapping was often carried out for individual cities by using high or medium spatial resolution images [12]. However, it is difficult to map urban areas using high or medium spatial resolution images such as Landsat TM/ETM+ at regional or global scales. First, it requires a large amount of human and computational resources to process and interpret these images. Second, obtaining a comprehensive set of good-quality images within any given year that covers a large study area is also difficult because of frequent cloudy conditions and infrequent measurements [14]. Therefore, developing new approaches for mapping urban land over large areas using coarse resolution images is an urgent task.

The Defense Meteorological Satellite Program’s Operational Line-scan System (DMSP-OLS) nighttime light (NTL) images have proven to be reliable sources of data for mapping urban areas at regional and global scales [10,15,16,17,18,19,20,21,22]. Threshold-based approaches are the most common techniques for extracting urban areas from DMSP-OLS data because of their simplicity [10,14,23,24]. Unfortunately, previous studies indicated that the threshold-based technique posed a problem to retrieve urban areas accurately. These approaches not only omit a large number of small settlements with low light brightness but also overestimate urban extents in larger-scale cities due to the blooming effect [9,10,24,25]. The multiple thresholds method was then developed to map urban areas at regional scales by considering the development level of cities [10,26]. A disadvantage of this method is that there are no general guidelines available for selecting threshold values based on different levels of socioeconomic development [25].

Additionally, the well-documented issue of saturation in bright urban cores means that NTL is unable to accurately estimate Impervious Surface Area (ISA) or separate different types of land cover [14,27]. Thus, uncertainty cannot be avoided when extracting urban area from DMSP-OLS data at regional and global scales using these types of methods. The combination of multi-source remotely sensed data might provide more information than any individual source because each source has different characteristics. Previous studies have demonstrated that the vegetation index was closely and inversely correlated with impervious surfaces [28,29]. Recently, a number of methods have attempted to combine the vegetation index and NTL data to obtain more accurate urban areas [14,21,25]. For example, Cao et al. (2009) [25] have shown the potential of combining DMSP-OLS and SPOT-VGT data for extracting urban extents using semi-automatic Support Vector Machines (SVM). However, the semi-automatic training process for SVM is both complex and time-consuming. Furthermore, small urban areas within mixed pixels, which are usually abundant in DMSP-OLS data, are also ignored. Hence, spatial information of towns and villages is inevitably omitted using the proposed SVM method. Lu et al. (2008) [14] proposed a Human Settlement Index (HSI) by integrating DMSP-OLS NTL with Normalized Difference Vegetation Index (NDVI) data derived from MODIS for mapping regional human settlements. However, HSI overcorrects the saturation in the peri-urban areas [21]. The vegetation adjusted NTL urban index (VANUI) was also developed to correct the NTL saturation [21]. However, NDVI is limited in its ability to distinguish urban lands from other non-urban types when those show the lowest values, such as water and barren lands. It has been proven that NDWI can effectively differentiate water bodies from non-water features [30]. NDWI can make up for the limitations of the HSI and the VANUI, which only consider the contribution of NDVI.

This paper proposes a novel index, called the Normalized Urban Areas Composite Index (NUACI), to obtain accurate and timely urban dynamics in large areas based on the combination of DMSP-OLS, the vegetation index, and the Normalized Difference Water Index (NDWI). The innovation of NUACI is that it can significantly reduce the blooming and saturation effects of NTL data in urban areas by incorporating NDWI, which improves its ability to distinguish urban lands from non-urban types with the lowest vegetation densities. China, an area experiencing unprecedented urbanization with rapid economic growth over the past three decades, was selected for the case study. Landsat TM/ETM+ images were used to extract urban areas as reference data. The DMSP-OLS, MODIS EVI and NDWI data were integrated to calculate NUACI values, which were then used to estimate ISA through regression models.

2. Pixel Blooming and Saturation for the DMSP-OLS Imagery

DMSP-OLS originated in the mid-1960s with the primary purpose of observing moonlit cloud cover [31]. It has low-light-sensing capabilities for detecting artificial lighting on the Earth’s surface at night and has been widely used in urban studies such as mapping urban areas, population estimation, carbon emissions inventory and economic activities, among others [14,32,33,34,35]. However, the well-documented issues of pixel blooming and saturation limit the ability for using NTL data to estimate urban areas accurately, which led to DMSP-OLS being relegated to more of a novelty than a scientific tool in the 1980s [16].

The spurious indication of light in a location without a light source is called pixel blooming, which is a widespread phenomenon in NTL data, especially in large metropolitan areas [24]. Three main causes for the pixel blooming phenomena in NTL images were summarized by Small et al. [24]. The first cause of pixel blooming is attributable to atmospheric scattering in the blue part of the spectrum (0.4 to 0.5 mm), which is used for low-light imaging on the DMSP-OLS. The coarse spatial resolution and scattered light result in blooming at the peripheries of light sources. Second, there is a large overlap in the footprints of adjacent pixels. Finally, the accumulation of geolocation errors in the compositing process is another factor that induces pixel blooming. The overlay images of Landsat TM/ETM+ and DMSP-OLS are shown in Figure 1 illustrate that a large number of unlighted areas are affected by the blooming effect. Agricultural lands, forests and water bodies are all depicted with light. Therefore, pixel blooming is an obstacle in using NTL images to estimate urban areas accurately. Many attempts have been made to reduce the blooming effect. Early effort is to use a low light threshold of detection frequency [16]. However, no single threshold is well suited to map urban areas for various cities [24]. Vegetation index is presented to reduce the overestimated spatial extent of lighted areas [14,21,25]. However, this method cannot distinguish water bodies near cities, such as coastal areas, from urban areas.

Another primary barrier to using NTL data for urban studies is pixel saturation, which refers to the fact that data values in urban core areas tend to be truncated because of the limited radiometric range of DMSP-OLS. The saturation phenomena are illustrated using digital number (DN) profiles across a latitudinal transect in Figure 2, which indicates that NTL values are invariant in urban areas due to saturation. In reality, the corrected NTL data without saturation would have some peaks in the urban core. Several techniques have been developed to correct saturation. For example, Ziskin et al. [36] attempted to calibrate radiance by varying the gain on the sensor. However, radiance calibrated nighttime light images were produced involving a significant amount of labor and cost. Recently, some simpler methods based on using the combination of MODIS NDVI and NTL data have been presented to mitigate the saturation effects, such as HSI and VANUI [14,21].

Figure 1. Overlay images of Landsat ETM+ and DMSP-OLS for the six selected cities.

Figure 1. Overlay images of Landsat ETM+ and DMSP-OLS for the six selected cities.

Figure 2. Latitudinal transects of calibrated NTL and saturated NTL in 2010 for Shenzhen, China.

Figure 2. Latitudinal transects of calibrated NTL and saturated NTL in 2010 for Shenzhen, China.

3. Normalized Urban Areas Composite Index (NUACI)

Previous studies have shown that vegetation abundance is highly negatively correlated with the distribution of impervious surfaces [29,37,38]. Hence, vegetation index can be used to assist in the estimation of urban areas and has also been applied to correct the blooming effect in HSI and VANUI indices. Actually, vegetation cover often shows discriminated patterns in different seasons, such as crops and grasslands [14]. Unfortunately, it is impractical to retrieve the distribution of urban areas using only a single-date vegetation index. EVI has proven to be insensitive to background noise. It avoids saturation at places with high densities of vegetation cover and more accurately reflects the true background condition of land cover distribution. The annual maximum enhanced vegetation index (EVI_max) derived from MODIS images can accurately provide information about vegetation covered areas and is certainly preferable for use as auxiliary data for mapping urban boundaries. However, the values of EVI_max for non-vegetation land cover, including water bodies, bare soils and barren lands, are similar to those of urban areas, making it difficult to distinguish urban areas from non-vegetation land using only the EVI_max index. NDWI can differentiate water bodies from non-water features [30]), so a combination of EVI_max and NDWI might be able to discriminate urban areas from vegetation land cover and water bodies. Therefore, EVI_max and NDWI are combined with NTL data in this paper to establish NUACI for extracting an urban fraction.

To validate the feasibility of integrating EVI_max and NDWI, 11 scenes of TM/ETM+ images in China were selected, from which we collected 4320 sample plots involving six types of land cover: water bodies, agricultural land, barren lands, forest, grass and urban. Figure 3 illustrates the correlation between EVI_max and NDWI for different land covers. The values of NDWI and EVI_max for different sample plots are generally scattered within a certain coordinate space. When considering only EVI_max for correcting the saturation and blooming effects in NTL images, we can find many sample plots—especially those belonging to urban lands and water bodies—that overlap along their vertical axes for the same value range of EVI_max. Urban lands cannot be easily separated from water bodies, barren lands or grasslands with low densities of vegetation cover. This phenomenon is usually typical for cities surrounded by barren lands or grass lands covered with low densities of vegetation, as illustrated in Figure 4. When NDWI is incorporated into the index for deriving urban information, different types of land cover show large differences between NDWI and EVI_max. Sample plots of different types are generally scattered within identifiable zones in the two-dimensional space constituted by NDWI and EVI_max. The sample plots belonging to urban lands are basically centralized within the circle regions having lower NDWI and EVI_max. Although the sample plots of water bodies present lower EVI_max, they generally show the highest NDWI, which makes them easy to separate from urban lands. In addition, barren lands with the lowest EVI_max have a lower NDWI than water bodies but a higher NDWI than other types, while farmlands have a higher EVI_max but lower NDWI. This characteristic pattern of differences makes urban lands easily distinguishable from farmlands, water bodies and barren lands, even though a few sample plots of barren lands and sparse grass fall into the circle region. Therefore, most non-urban areas with comparatively higher DN values in NTL images can be excluded from urban lands by incorporating both EVI_max and NDWI.

Figure 3. Relationship between NDWI and EVI_max for sample plots with different land cover types.

Figure 3. Relationship between NDWI and EVI_max for sample plots with different land cover types.

As described above, because different characteristics can be obtained from the DMSP-OLS, MODIS EVI_max and NDWI datasets, it is necessary to identify the relationships among those datasets to improve the accuracy of mapping urban areas. Therefore, a new index, the NUACI, is established in this paper, as expressed in Equation (1).

$$\text{NUACI}=\{\begin{array}{l}0, d>r, d=\sqrt{{(\text{NDWI}-a)}^2+{({\text{EVI}}_{\max }-b)}^2}\\ \left(1-d/r\right)\ast \frac{{\text{OLS-OLS}}_{\min }}{{\text{OLS}}_{\max }-{\text{OLS}}_{\min }} , d\le r\end{array}$$

(1)

where a and b are the average value of NDWI and EVI_max for urban samples, respectively, r is the radius of the circle region aggregated with urban samples, d is the distance to the center of the circle and OLS_min and OLS_max are the minimum and maximum values in the DMSP-OLS image, respectively. The value of NUACI can be standardized within a range varying from 0 to 1. The true distribution pattern of urban lands can be reflected according to the adoption value of r. NDWI can be calculated as follows:

$$\text{NDWI}=\frac{{\rho }_{857}-{\rho }_{1241}}{{\rho }_{857}+{\rho }_{1241}}$$

(2)

where ρ₈₅₇ is reflectance in the near infrared band at 857 nm, and ρ₁₂₄₁ is reflectance in the water absorption band at 1241 nm. EVI can be expressed as follows:

$$\text{EVI}=\frac{{\rho }_N^{*}-{\rho }_R^{*}}{{\rho }_N^{*}+C_1{\rho }_R^{*}-C_2{\rho }_B^{*}+L}$$

(3)

where ${\rho }_N^{*}$ is reflectance in the near infrared band, ${\rho }_R^{*}$ is reflectance in the red band, L is the coefficient of soil adaption and C₁ and C₂ are the coefficients used to correct the influence of atmospheric aerosol scattering in the blue and red bands, respectively.

Figure 4. A typical spatial pattern of land cover (e.g., grass covered with low density of vegetation around urban land) in a city.

Figure 4. A typical spatial pattern of land cover (e.g., grass covered with low density of vegetation around urban land) in a city.

4. Materials and Data Preprocessing

The whole of China was selected for the case study. The country has experienced unprecedented urbanization since the 1978 reforms. A time series of DMSP-OLS NTL images, Terra MODIS reflectance and EVI products as well as Landsat ETM+ were used in this study. DMSP-OLS version 4 NTLs datasets with 1 km spatial resolution from 2000, 2005 and 2010 were downloaded from the National Oceanic and Atmospheric Administration (NOAA)/National Geophysical Data Center (NGDC) website. Each dataset is composed of average DN values of stable light images for an entire year. DN values for lights on NTL images range from 1 to 63, while the background and noise are all recorded as zero. The subsets of NTL images were then extracted according to Chinese administrative boundaries.

Surface reflectance products of Terra MODIS (MOD09A1) and EVI (MOD13A1) with 463 m spatial resolution were downloaded from the National Aeronautics and Space Administration (NASA) website. MOD09A1 products, which were preprocessed to include atmospheric, radiometric and geometric correction, were used to calculate the NDWI. The 16-day EVI product was processed by compositing maximum values to create a yearly maximal EVI image (EVI_max). The DMSP-OLS, NDWI and EVI_max images were co-registered and projected to the Albers Conical Equal Area projection. All the datasets were resampled to the same spatial resolution of 500 m.

Thirty scenes of Landsat ETM+ images covering partial regions of China with spatial resolutions of 28.5 m from the years 2000, 2005 and 2010 were also obtained from the United States Geological Survey website and used to extract actual impervious surface area, which we then used as reference data to establish regression models and validate the performance of NUACI. All the selected TM/ETM+ images were reprojected to the Albers Conical Equal Area projection. The urban areas were extracted using the maximum likelihood classifier (MLC), which has been commonly applied to the classification of medium resolution remote sensing images. Table 1 shows the confusion matrices of the urban land classifications extracted from Landsat ETM+ images. In Table 1, the overall accuracy shows high values of 91.13% (2000), 91.58% (2005), and 91.55% (2010), which indicate that the classification image is fairly accurate. Therefore, the urban information derived from the TM/ETM+ images can be used as reliable reference data. The obtained urban classification was then aggregated to calculate actual ISA with a pixel size of 500 m by 500 m, which matches the spatial resolution of MODIS EVI, NDWI, and DMSP-OLS data.

Table 1. Confusion matrices for urban land classifications (2000, 2005, 2010).

Classification	Reference
	2000		2005		2010
	Urban	Non-urban	Urban	Non-urban	Urban	Non-urban
Urban	1141	163	1186	153	1239	169
Nonurban	447	5123	426	5109	412	5054
Over Accuracy	91.13%		91.58%		91.55%

Table 1. Confusion matrices for urban land classifications (2000, 2005, 2010).

Classification	Reference
	2000		2005		2010
	Urban	Non-urban	Urban	Non-urban	Urban	Non-urban
Urban	1141	163	1186	153	1239	169
Nonurban	447	5123	426	5109	412	5054
Over Accuracy	91.13%		91.58%		91.55%

5. Model Application and Results Discussion

5.1. Calculation and Validation of NUACI

Nineteen MODIS surface reflectance data and 16-day composite EVI products covering China were used to obtain NDWI and EVI_max, respectively. The NUACI for the whole of China was then calculated using the combination of DMSP-OLS, NDWI and EVImax images with ArcMap model builder, which is used to implement GIS processes by linking input data automatically. The parameters a, b and r in Equation (1) must be determined before implementing NUACI calculation for the study region. Urban samples collected from ETM+ can be used to determine the a, b and r parameters. We calculated the average values of NDWI and EVI_max for urban samples to obtain the values of a and b. The parameter r can be set by calculating the distance between the farthest urban sample and the central urban sample, the position of which is given as (a, b). In this paper, the parameters a, b and r were finally set to −0.35, 0.15 and 0.4, respectively. Six cities undergoing rapid urbanization, including Beijing, Tianjin, Hangzhou, Guangzhou, Shenzhen, and Wuhan, were selected to make the comparison between the NUACI data and the NTL images. As shown in Figure 5, NTL values were kept invariant in the urban core because of pixel saturation. Furthermore, the blooming effect caused pixels scattered near urban periphery to exhibit high DN values in NTL images. Non-urban land cover within peri-urban areas can be effectively separated from urban lands by using NUACI.

To further evaluate the ability of NUACI on reducing pixel saturation and eliminating the blooming effect, a visual comparison was performed between NUACI, NTL, calibrated NTL without saturation, VANUI and HSI for a latitudinal transect in Shenzhen (Figure 6). It was apparent that the spatial pattern of NUACI is similar to those of VANUI and HSI. However, NUACI values vary within the urban core, increase toward the urban core, and decrease away from the urban core. The comparison illustrated that although these three indices can reduce pixel saturation and increase inter-urban variability, NUACI values have more variations within urban areas. More importantly, both VANUI and HSI are incapable of reducing the blooming effect because their values for non-urban land cover around the peri-urban areas are still high. In contrast, NUACI values for those areas are very low—almost equal to zero. This result indicates that NUACI has the ability to solve the pixel blooming problem. Figure 7 and Figure 8 also display visual comparisons among corrected images selected from the Pearl Delta and from the city of Wuhan as sample areas. The shoreline or riverbank located in the cities is apparently affected by pixel blooming in the NTL data. Pixel blooming is most evident on the imagery where glints of light into adjacent water bodies can be observed. As illustrated in Figure 7 and Figure 8, both VANUI and HSI are incapable of separating water bodies from urban areas, whereas NUACI can accurately capture the spatial distribution information of urban areas.

Figure 5. Comparison of Landsat ETM+ images (a); NTL data (b); and NUACI images (c) for six cities.

Figure 5. Comparison of Landsat ETM+ images (a); NTL data (b); and NUACI images (c) for six cities.

Further quantitative assessment was carried out to validate the performance of the proposed NUACI index on eliminating the blooming effect. In this experiment, the Pearl River Delta was selected as the sample area. The values of areas located within present urban lands were discussed and analyzed in detail. Urban lands were retrieved from TM/ETM+ images covering the Pearl River Delta, a total of six buffers around existing urban areas within an interval distance of 0.5 km were designated. The mean values of those buffer areas were then calculated for different urban indices. Theoretically, DN values for those non-urban areas should equal zero. However, the buffer areas may show different mean values because of the blooming effect. When lower mean values appear in the corrected images, we can consider the blooming effect extremely restrained. Table 2 shows the average values for different indices based on images within six buffers. The lowest average values were obtained from the NUACI-based images, which were less than those from DMSP-OLS, HSI, and VANUI images in all the buffer areas.

Figure 9 shows a visual comparison of value distribution pattern among the different urban indices images. The scope of blooming effects is obviously reduced and sometimes even eliminated in the NUACI images. When the buffer area contains large-scale water bodies, the values of NUACI are approximately zero, whereas pixels within the buffer areas for the other urban indices images show high values because of the remaining blooming effect.

Figure 6. Latitudinal transects of calibrated NTL, saturated NTL, NUACI, HSI, and VANUI for Shenzhen, China.

Figure 6. Latitudinal transects of calibrated NTL, saturated NTL, NUACI, HSI, and VANUI for Shenzhen, China.

Figure 7. Comparison of different indexes. (a) Landsat ETM+ images; (b) saturated NTL; (c) calibrated NTL; (d) VANUI; (e) HIS; and (f) NUACI for the Pearl River Delta, China.

Figure 7. Comparison of different indexes. (a) Landsat ETM+ images; (b) saturated NTL; (c) calibrated NTL; (d) VANUI; (e) HIS; and (f) NUACI for the Pearl River Delta, China.

Figure 8. Comparison of different indexes. (a) Landsat ETM+ images; (b) saturated NTL; (c) calibrated NTL; (d) VANUI; (e) HIS; and (f) NUACI for Wuhan, China.

Figure 8. Comparison of different indexes. (a) Landsat ETM+ images; (b) saturated NTL; (c) calibrated NTL; (d) VANUI; (e) HIS; and (f) NUACI for Wuhan, China.

Figure 9. Comparison of Landsat TM/ETM+ images, saturated NTL, VANUI, HSI, and NUACI for buffer areas: (a) Dongguan; (b) Zhuhai; (c) Guangzhou; (d) Shenzhen; and (e) Huizhou.

Figure 9. Comparison of Landsat TM/ETM+ images, saturated NTL, VANUI, HSI, and NUACI for buffer areas: (a) Dongguan; (b) Zhuhai; (c) Guangzhou; (d) Shenzhen; and (e) Huizhou.

Table 2. Comparison of the average values within different buffer areas.

City	0.5 km				1 km				1.5 km
	OLS	HSI	VANUI	NUACI	OLS	HSI	VANUI	NUACI	OLS	HSI	VANUI	NUACI
Foshan	0.5326	1.4353	0.3640	0.2669	0.4771	1.2940	0.3236	0.2297	0.4424	1.2120	0.2987	0.2087
Guangzhou	0.4042	0.9953	0.2472	0.1364	0.3442	0.8680	0.2069	0.1037	0.3101	0.8082	0.1857	0.0886
Huizhou	0.2257	0.5773	0.1220	0.0431	0.1816	0.5108	0.0957	0.0290	0.1541	0.4766	0.0804	0.0223
Dongguan	0.8232	1.9550	0.5233	0.3224	0.7883	1.8401	0.4934	0.2799	0.7597	1.7724	0.4730	0.2566
Shenzhen	0.7572	1.9056	0.4385	0.1823	0.7064	2.0489	0.4104	0.1452	0.6696	2.2156	0.3977	0.1294
Zhongshan	0.6499	1.3683	0.3945	0.2301	0.6060	1.2739	0.3651	0.2042	0.5767	1.2167	0.3460	0.1889
Zhuhai	0.3360	0.9272	0.2052	0.0870	0.3030	0.9436	0.1903	0.0742	0.2777	0.9643	0.1791	0.0659
City	2 km				2.5 km				5 km
	OLS	HSI	VANUI	NUACI	OLS	HSI	VANUI	NUACI	OLS	HSI	VANUI	NUACI
Foshan	0.4229	1.1671	0.2846	0.1977	0.4083	1.1315	0.2739	0.1894	0.3571	1.0127	0.2371	0.1610
Guangzhou	0.2865	0.7673	0.1713	0.0799	0.2693	0.7353	0.1606	0.0743	0.2169	0.6367	0.1278	0.0578
Huizhou	0.1360	0.4563	0.0707	0.0188	0.1220	0.4408	0.0633	0.0163	0.0886	0.4063	0.0457	0.0109
Dongguan	0.7360	1.7257	0.4578	0.2423	0.7154	1.6886	0.4453	0.2321	0.6624	1.5885	0.4136	0.2091
Shenzhen	0.6427	2.3297	0.3894	0.1166	0.6242	2.4109	0.3857	0.1083	0.5580	2.3948	0.3629	0.0877
Zhongshan	0.5543	1.1837	0.3326	0.1792	0.5362	1.1639	0.3224	0.1722	0.4904	1.1131	0.2958	0.1521
Zhuhai	0.2548	0.9892	0.1681	0.0586	0.2360	1.0119	0.1586	0.0532	0.1669	1.0565	0.1175	0.0333

Table 2. Comparison of the average values within different buffer areas.

City	0.5 km				1 km				1.5 km
	OLS	HSI	VANUI	NUACI	OLS	HSI	VANUI	NUACI	OLS	HSI	VANUI	NUACI
Foshan	0.5326	1.4353	0.3640	0.2669	0.4771	1.2940	0.3236	0.2297	0.4424	1.2120	0.2987	0.2087
Guangzhou	0.4042	0.9953	0.2472	0.1364	0.3442	0.8680	0.2069	0.1037	0.3101	0.8082	0.1857	0.0886
Huizhou	0.2257	0.5773	0.1220	0.0431	0.1816	0.5108	0.0957	0.0290	0.1541	0.4766	0.0804	0.0223
Dongguan	0.8232	1.9550	0.5233	0.3224	0.7883	1.8401	0.4934	0.2799	0.7597	1.7724	0.4730	0.2566
Shenzhen	0.7572	1.9056	0.4385	0.1823	0.7064	2.0489	0.4104	0.1452	0.6696	2.2156	0.3977	0.1294
Zhongshan	0.6499	1.3683	0.3945	0.2301	0.6060	1.2739	0.3651	0.2042	0.5767	1.2167	0.3460	0.1889
Zhuhai	0.3360	0.9272	0.2052	0.0870	0.3030	0.9436	0.1903	0.0742	0.2777	0.9643	0.1791	0.0659
City	2 km				2.5 km				5 km
	OLS	HSI	VANUI	NUACI	OLS	HSI	VANUI	NUACI	OLS	HSI	VANUI	NUACI
Foshan	0.4229	1.1671	0.2846	0.1977	0.4083	1.1315	0.2739	0.1894	0.3571	1.0127	0.2371	0.1610
Guangzhou	0.2865	0.7673	0.1713	0.0799	0.2693	0.7353	0.1606	0.0743	0.2169	0.6367	0.1278	0.0578
Huizhou	0.1360	0.4563	0.0707	0.0188	0.1220	0.4408	0.0633	0.0163	0.0886	0.4063	0.0457	0.0109
Dongguan	0.7360	1.7257	0.4578	0.2423	0.7154	1.6886	0.4453	0.2321	0.6624	1.5885	0.4136	0.2091
Shenzhen	0.6427	2.3297	0.3894	0.1166	0.6242	2.4109	0.3857	0.1083	0.5580	2.3948	0.3629	0.0877
Zhongshan	0.5543	1.1837	0.3326	0.1792	0.5362	1.1639	0.3224	0.1722	0.4904	1.1131	0.2958	0.1521
Zhuhai	0.2548	0.9892	0.1681	0.0586	0.2360	1.0119	0.1586	0.0532	0.1669	1.0565	0.1175	0.0333

Additionally, urban indices such as HSI and VANUI have been directly applied to extract the extent of urban areas by using the threshold segmentation method. The accuracy of urban land extraction can reflect the efficiency on reducing the blooming effect. In this paper, we extracted urban areas using the proposed NUACI. We selected a total of eleven cities as samples. Urban areas obtained from normalized DMSP-OLS, HSI, and VANUI images were further evaluated to validate the efficiency of NUACI in eliminating the blooming effect. In addition, urban land extracted from the MODIS land-use products with a resolution of 500 m, which were generated using SVM methods [2], were also downloaded for comparison purposes. The total area of urban lands was derived from Landsat TM/ETM+ images using the MLC method and used as actual reference data. Corresponding equal areas of urban land for each sample city were retrieved from DMSP-OLS NTL data, MODIS land-use data, HSI image, VANUI image, and NUACI image, respectively. Quantitative indices including Kappa coefficient and overall accuracy were calculated by means of point-to-point comparison. The results were applied to measure the performance of deriving urban areas with these different combinations.

Figure 10, Figure 11 and Figure 12 show the spatial distribution of classification accuracy, commission errors, and omission errors on urban areas derived from the different images from 2000, 2005, and 2010. All the results derived from NUACI exhibit lower commission and omission errors in all the sample areas. The spatial patterns are generally similar to those derived from MODIS land-use data. Within the time variation, the spatial patterns of urban lands obtained with NUACI method for the selected cities are more rational than those from MODIS land-use data. Higher commission errors (in green) generally presented in the urban areas extracted from DMSP-OLS NTL because of the blooming effect on the coastal area, rivers, vegetation covered lands, barren lands and water bodies. The results derived from HSI and VANUI can reduce the blooming effect on vegetation covered areas, which also reduces commission errors on vegetation covered lands. In addition, the blooming effect on coastal areas, rivers and large-scale water bodies are further efficiently reduced by using NUACI, although a few commission errors were distributed along small-scale rivers. The densities of omission pixels also decreased in the extracted urban areas when using the NUACI method. In addition, results derived from MODIS land-use data presented higher omission errors than those from NUACI in regions with low urban land densities. As shown in Table 3, Table 4 and Table 5, almost all the cities show the highest Kappa coefficients and overall accuracy of urban areas when extracted from NUACI compared with those from DMSP-OLS, HSI, VANUI, and the MODIS data.

Figure 10. Comparison of the results of extracting urban lands by the Threshold Segmentation Method in the year 2000: (a) Guangzhou; (b) Shanghai; (c) Shenzhen; and (d) Hangzhou.

Figure 10. Comparison of the results of extracting urban lands by the Threshold Segmentation Method in the year 2000: (a) Guangzhou; (b) Shanghai; (c) Shenzhen; and (d) Hangzhou.

Figure 11. Comparison of the results of extracting urban lands by the Threshold Segmentation Method in the year 2005: (a) Guangzhou; (b) Shanghai; (c) Shenzhen; and (d) Hangzhou.

Figure 11. Comparison of the results of extracting urban lands by the Threshold Segmentation Method in the year 2005: (a) Guangzhou; (b) Shanghai; (c) Shenzhen; and (d) Hangzhou.

Figure 12. Comparison of the results of extracting urban lands by the Threshold Segmentation Method in the year 2010: (a) Guangzhou; (b) Shanghai; (c) Shenzhen; and (d) Hangzhou.

Figure 12. Comparison of the results of extracting urban lands by the Threshold Segmentation Method in the year 2010: (a) Guangzhou; (b) Shanghai; (c) Shenzhen; and (d) Hangzhou.

Table 3. Kappa and overall accuracy of extracting urban land in 2000.

Name	Urban Area (KM2)	Kappa					Overall Accuracy
		DMSP	HSI	VANUI	MODIS 500 m	NUACI	DMSP	HSI	VANUI	MODIS 500 m	NUACI
Beijing	1047	0.5550	0.6329	0.6224	0.5567	0.6393	0.7795	0.8181	0.8129	0.7817	0.8213
Shanghai	760	0.5432	0.5376	0.4496	0.5549	0.6279	0.7720	0.7692	0.7253	0.7801	0.8143
Guangzhou	546	0.4935	0.4667	0.3786	0.4459	0.4916	0.8013	0.7851	0.7384	0.7517	0.8023
Shenzhen	606	0.4404	0.4763	0.2946	0.5245	0.5931	0.7367	0.7682	0.6709	0.7966	0.8344
Wuhan	323	0.5078	0.5189	0.3564	0.6202	0.6009	0.8383	0.8419	0.7488	0.8636	0.8689
Nanjing	257	0.4660	0.4778	0.3923	0.4369	0.5566	0.7396	0.7454	0.7037	0.7399	0.7838
Chongqing	204	0.5417	0.5701	0.5193	0.5188	0.5739	0.8396	0.8495	0.8318	0.8262	0.8508
Dalian	281	0.5712	0.4651	0.0081	0.6293	0.6939	0.8858	0.8491	0.4029	0.8986	0.9186
Hangzhou	174	0.3606	0.3728	0.3587	0.5104	0.4548	0.7427	0.7481	0.7423	0.8044	0.7808
Changsha	129	0.4229	0.4369	0.3782	0.4677	0.4398	0.7714	0.7769	0.7537	0.7786	0.7781
Xiamen	107	0.2808	0.1655	−0.0090	0.3012	0.4070	0.8370	0.6967	0.3913	0.7669	0.8660

Table 3. Kappa and overall accuracy of extracting urban land in 2000.

Name	Urban Area (KM2)	Kappa					Overall Accuracy
		DMSP	HSI	VANUI	MODIS 500 m	NUACI	DMSP	HSI	VANUI	MODIS 500 m	NUACI
Beijing	1047	0.5550	0.6329	0.6224	0.5567	0.6393	0.7795	0.8181	0.8129	0.7817	0.8213
Shanghai	760	0.5432	0.5376	0.4496	0.5549	0.6279	0.7720	0.7692	0.7253	0.7801	0.8143
Guangzhou	546	0.4935	0.4667	0.3786	0.4459	0.4916	0.8013	0.7851	0.7384	0.7517	0.8023
Shenzhen	606	0.4404	0.4763	0.2946	0.5245	0.5931	0.7367	0.7682	0.6709	0.7966	0.8344
Wuhan	323	0.5078	0.5189	0.3564	0.6202	0.6009	0.8383	0.8419	0.7488	0.8636	0.8689
Nanjing	257	0.4660	0.4778	0.3923	0.4369	0.5566	0.7396	0.7454	0.7037	0.7399	0.7838
Chongqing	204	0.5417	0.5701	0.5193	0.5188	0.5739	0.8396	0.8495	0.8318	0.8262	0.8508
Dalian	281	0.5712	0.4651	0.0081	0.6293	0.6939	0.8858	0.8491	0.4029	0.8986	0.9186
Hangzhou	174	0.3606	0.3728	0.3587	0.5104	0.4548	0.7427	0.7481	0.7423	0.8044	0.7808
Changsha	129	0.4229	0.4369	0.3782	0.4677	0.4398	0.7714	0.7769	0.7537	0.7786	0.7781
Xiamen	107	0.2808	0.1655	−0.0090	0.3012	0.4070	0.8370	0.6967	0.3913	0.7669	0.8660

Table 4. Kappa and overall accuracy of extracting urban land in 2005.

Name	Urban Area (km2)	Kappa					Overall Accuracy
		DMSP	HSI	VANUI	MODIS 500 m	NUACI	DMSP	HSI	VANUI	MODIS 500 m	NUACI
Beijing	1236	0.5610	0.5989	0.5994	0.5020	0.6058	0.7993	0.8168	0.8168	0.7625	0.8198
Shanghai	882	0.4816	0.4564	0.5103	0.5909	0.6031	0.7523	0.7399	0.7657	0.8080	0.8101
Guangzhou	799	0.4566	0.5134	0.5244	0.4934	0.5420	0.7421	0.7691	0.7743	0.7593	0.7828
Shenzhen	813	0.4150	0.4276	0.4841	0.4694	0.5979	0.7238	0.7300	0.7566	0.7559	0.8103
Wuhan	417	0.5030	0.5107	0.5271	0.5991	0.5791	0.8052	0.8083	0.8147	0.8444	0.8351
Nanjing	337	0.4501	0.4192	0.4532	0.2861	0.5001	0.7282	0.7130	0.7298	0.6232	0.7529
Chongqing	250	0.4817	0.5395	0.5342	0.4940	0.5609	0.7921	0.8151	0.8132	0.8042	0.8238
Dalian	296	0.5694	0.3849	0.5095	0.6290	0.6534	0.8804	0.8293	0.8639	0.8964	0.9039
Hangzhou	238	0.3369	0.3915	0.4003	0.4665	0.4903	0.6871	0.7127	0.7169	0.7618	0.7594
Changsha	170	0.4165	0.4406	0.4410	0.4317	0.4570	0.7318	0.7428	0.7430	0.7458	0.7503
Xiamen	166	0.3261	0.0188	0.1469	0.4007	0.4841	0.7835	0.6851	0.7260	0.7810	0.8345

Table 4. Kappa and overall accuracy of extracting urban land in 2005.

Name	Urban Area (km2)	Kappa					Overall Accuracy
		DMSP	HSI	VANUI	MODIS 500 m	NUACI	DMSP	HSI	VANUI	MODIS 500 m	NUACI
Beijing	1236	0.5610	0.5989	0.5994	0.5020	0.6058	0.7993	0.8168	0.8168	0.7625	0.8198
Shanghai	882	0.4816	0.4564	0.5103	0.5909	0.6031	0.7523	0.7399	0.7657	0.8080	0.8101
Guangzhou	799	0.4566	0.5134	0.5244	0.4934	0.5420	0.7421	0.7691	0.7743	0.7593	0.7828
Shenzhen	813	0.4150	0.4276	0.4841	0.4694	0.5979	0.7238	0.7300	0.7566	0.7559	0.8103
Wuhan	417	0.5030	0.5107	0.5271	0.5991	0.5791	0.8052	0.8083	0.8147	0.8444	0.8351
Nanjing	337	0.4501	0.4192	0.4532	0.2861	0.5001	0.7282	0.7130	0.7298	0.6232	0.7529
Chongqing	250	0.4817	0.5395	0.5342	0.4940	0.5609	0.7921	0.8151	0.8132	0.8042	0.8238
Dalian	296	0.5694	0.3849	0.5095	0.6290	0.6534	0.8804	0.8293	0.8639	0.8964	0.9039
Hangzhou	238	0.3369	0.3915	0.4003	0.4665	0.4903	0.6871	0.7127	0.7169	0.7618	0.7594
Changsha	170	0.4165	0.4406	0.4410	0.4317	0.4570	0.7318	0.7428	0.7430	0.7458	0.7503
Xiamen	166	0.3261	0.0188	0.1469	0.4007	0.4841	0.7835	0.6851	0.7260	0.7810	0.8345

Table 5. Kappa and overall accuracy of extracting urban land in 2010.

Name	Urban Area (km2)	Kappa					Overall Accuracy
		DMSP	HSI	VANUI	MODIS 500 m	NUACI	DMSP	HSI	VANUI	MODIS 500 m	NUACI
Beijing	1675	0.4718	0.4587	0.3900	0.1944	0.4861	0.8852	0.8824	0.8674	0.6448	0.8717
Shanghai	1017	0.4460	0.3454	0.2856	0.5448	0.5836	0.7946	0.7232	0.6980	0.7989	0.8241
Guangzhou	925	0.4530	0.4140	0.2879	0.4741	0.4820	0.7296	0.7103	0.6479	0.7432	0.7438
Shenzhen	969	0.4225	0.4319	0.2888	0.4690	0.6584	0.7102	0.7183	0.6473	0.7427	0.8306
Wuhan	581	0.4454	0.3918	0.2632	0.5454	0.4976	0.7407	0.7157	0.6468	0.8006	0.7651
Nanjing	407	0.3255	0.4240	0.3669	0.1918	0.5651	0.7011	0.7448	0.7195	0.5273	0.8073
Chongqing	335	0.4649	0.5314	0.5450	0.4354	0.5560	0.7501	0.7811	0.7875	0.7548	0.7926
Dalian	368	0.6169	0.4505	0.0366	0.6203	0.6914	0.8658	0.7843	0.4127	0.8837	0.8987
Hangzhou	256	0.3362	0.3904	0.3874	0.4446	0.5061	0.6778	0.7048	0.7033	0.7446	0.7608
Changsha	263	0.4335	0.4430	0.4026	0.3281	0.4885	0.7202	0.7250	0.7049	0.6511	0.7474
Xiamen	231	0.3015	0.0730	-0.0636	0.4348	0.5964	0.7182	0.5220	0.3787	0.7735	0.8373

Table 5. Kappa and overall accuracy of extracting urban land in 2010.

Name	Urban Area (km2)	Kappa					Overall Accuracy
		DMSP	HSI	VANUI	MODIS 500 m	NUACI	DMSP	HSI	VANUI	MODIS 500 m	NUACI
Beijing	1675	0.4718	0.4587	0.3900	0.1944	0.4861	0.8852	0.8824	0.8674	0.6448	0.8717
Shanghai	1017	0.4460	0.3454	0.2856	0.5448	0.5836	0.7946	0.7232	0.6980	0.7989	0.8241
Guangzhou	925	0.4530	0.4140	0.2879	0.4741	0.4820	0.7296	0.7103	0.6479	0.7432	0.7438
Shenzhen	969	0.4225	0.4319	0.2888	0.4690	0.6584	0.7102	0.7183	0.6473	0.7427	0.8306
Wuhan	581	0.4454	0.3918	0.2632	0.5454	0.4976	0.7407	0.7157	0.6468	0.8006	0.7651
Nanjing	407	0.3255	0.4240	0.3669	0.1918	0.5651	0.7011	0.7448	0.7195	0.5273	0.8073
Chongqing	335	0.4649	0.5314	0.5450	0.4354	0.5560	0.7501	0.7811	0.7875	0.7548	0.7926
Dalian	368	0.6169	0.4505	0.0366	0.6203	0.6914	0.8658	0.7843	0.4127	0.8837	0.8987
Hangzhou	256	0.3362	0.3904	0.3874	0.4446	0.5061	0.6778	0.7048	0.7033	0.7446	0.7608
Changsha	263	0.4335	0.4430	0.4026	0.3281	0.4885	0.7202	0.7250	0.7049	0.6511	0.7474
Xiamen	231	0.3015	0.0730	-0.0636	0.4348	0.5964	0.7182	0.5220	0.3787	0.7735	0.8373

5.2. Estimation of ISA with NUACI

The quantitative and visual comparison validated incorporating EVI and NDWI to efficiently reduce the blooming effect and saturation limitation. NUACI has the ability to retrieve accurate urban information using the coarse resolution NTL data. In actuality, mapping the national urban fraction is more attractive and has increased interest in using the coarse resolution images. ISA has been considered as an important indicator for measuring the degree of urbanization. Numerous studies have shown that NTL data were positively correlated with ISA data [14,21,39]. A regression model was thus established with the NUACI for estimating the density of ISA for the whole of China. In the linear regression model, the aggregated ISA data derived from the retrieved thirty scenes of Landsat TM/ETM+ classification results were used as reference data (dependent variable), while NUACI was used as a single independent variable. ISA extracted from regression models based on normalized DMSP-OLS, VANUI and HSI were used in the comparison to evaluate the ability of NUACI to extract urban fractions. Correlation coefficients between the estimated ISA and the reference data were used to examine the performance of the regression models. In addition, the root-mean square error (RMSE), which is a quadratic scoring rule available to measure the average magnitude of the error, was used to perform an accuracy assessment. RMSE can be calculated as follows:

$$RMSE=\sqrt{(1/N){\displaystyle \sum _{i=1}^N}{(p_i-t_i)}^2}$$

(4)

where N is the number of samples, and p_i and t_i are the measured and estimated values, respectively.

Five thousand sample plots were randomly selected from the reference data. From those, 3000 samples were used as the training data to establish regression models corresponding to normalized DMSP-OLS, NUACI, VANUI and HSI, respectively. The remaining 2000 samples were used to calculate quantitative indices including the coefficient of determination (R²), the correlation coefficient (R) and RMSE for validating the regression models. The estimation and validation were performed using Matlab, a programming language developed by MathWorks [40]. Table 6 lists the regression models and the calculated evaluation indices. It shows that the R² value of regression models based on DMSP-OLS, HSI, VANUI, and NUACI are 0.4459, 0.6888, 0.7822, and 0.8079, respectively. The regression models using VANUI or NUACI as the independent variable provide comparatively higher R² values than those based on normalized DMSP-OLS or HSI. The NUACI-based regression model has the highest R² value of 0.8079. This result implies that NUACI could be the most suitable method for estimating ISA. In contrast, the DMSP-OLS-based regression model has the lowest R² value of 0.4459. The results for correlation coefficients (R) were similar to the R² values. The NUACI-based regression model has the highest R value (0.8908), while the DMSP-OLS, HSI and VANUI based regression models had R values of 0.6680, 0.8241 and 0.8695, respectively. This implies that the NUACI-based regression model has the best performance for testing samples. Moreover, the RMSE of the NUACI-based regression model shows the lowest value of 0.1176. The scatterplots (Figure 13) generated using the estimations and the reference data also demonstrate that there is higher correlation between ISA and NUACI or VANUI than between ISA and other urban indices. In addition, the three regression models other than the NUACI model are inclined to overestimate the fraction values for some pixels that have zero value (horizontal axis). This indicates that the NUACI-based regression model can successfully reduce the blooming effect and achieve better estimation accuracy of ISA than other models.

Table 6. Regression models based on different variables for estimating ISA.

Variable Used	Regression Model (y = ax + b)	R2	R	RMSE
DMSP	ISA = 2.589DMSP − 0.159	0.4459	0.6680	0.1667
HSI	ISA = 0.669HSI − 0.825	0.6888	0.8241	0.1541
VANUI	ISA = 3.222VANUI − 1.307	0.7822	0.8695	0.1458
NUACI	ISA = 1.701NUACI − 0.261	0.8079	0.8908	0.1176

Table 6. Regression models based on different variables for estimating ISA.

Variable Used	Regression Model (y = ax + b)	R2	R	RMSE
DMSP	ISA = 2.589DMSP − 0.159	0.4459	0.6680	0.1667
HSI	ISA = 0.669HSI − 0.825	0.6888	0.8241	0.1541
VANUI	ISA = 3.222VANUI − 1.307	0.7822	0.8695	0.1458
NUACI	ISA = 1.701NUACI − 0.261	0.8079	0.8908	0.1176

According to the comparison above, NUACI can perform efficiently for retrieving the ISA distribution in the sample area. Therefore, the NUACI-based regression models were applied to the estimation of ISA in 2000, 2005 and 2010 for the whole of China. As shown in Table 7, there are some variations in the parameters of regression models for different years. The regression model in 2005 has the highest estimation accuracy. Overall, a large proportion of ISA are distributed in the eastern regions, especially in the Pearl River Delta, Yangtze River Delta and Beijing-Tianjin-Tangshan metropolitan areas as shown in Figure 14. The ISA growth distributions from 2000 to 2010 are shown in Figure 15, which illustrates that the Pearl River Delta, Yangtze River Delta and Beijing-Tianjin-Tangshan metropolitan areas are the three major regions with high ISA growth rates. In particular, the highest ISA expansion rates occur in the Yangtze River Delta. Figure 15 also demonstrates that medium and small cities surrounding megacities have rapid ISA growth rates because of the radiation effect.

Figure 13. Comparison of the scatterplots between the estimated ISA and the reference data.

Figure 13. Comparison of the scatterplots between the estimated ISA and the reference data.

Table 7. Regression models for estimating ISA based on NUACI in 2000, 2005 and 2010.

Year	Regression Model (y = ax + b)	R2	R	RMSE
2000	ISA = 1.750*NUACI − 0.430	0.7914	0.8413	0.1400
2005	ISA = 1.701*NUACI − 0.261	0.8079	0.8908	0.1176
2010	ISA = 1.669*NUACI − 0.121	0.7834	0.7382	0.1393

Table 7. Regression models for estimating ISA based on NUACI in 2000, 2005 and 2010.

Year	Regression Model (y = ax + b)	R2	R	RMSE
2000	ISA = 1.750*NUACI − 0.430	0.7914	0.8413	0.1400
2005	ISA = 1.701*NUACI − 0.261	0.8079	0.8908	0.1176
2010	ISA = 1.669*NUACI − 0.121	0.7834	0.7382	0.1393

Figure 14. Estimated ISA of the three metropolitan areas in the years 2000, 2005 and 2010 using the NUACI-based regression model: (a–c) ISA estimated for the Beijing-Tianjin-Tangshan metropolitan areas in 2000, 2005, and 2010; (d–f) ISA estimated for the Yangtze River Delta in 2000, 2005, and 2010; and (g–i) ISA estimated for the Pearl River Delta in 2000, 2005, and 2010.

Figure 14. Estimated ISA of the three metropolitan areas in the years 2000, 2005 and 2010 using the NUACI-based regression model: (a–c) ISA estimated for the Beijing-Tianjin-Tangshan metropolitan areas in 2000, 2005, and 2010; (d–f) ISA estimated for the Yangtze River Delta in 2000, 2005, and 2010; and (g–i) ISA estimated for the Pearl River Delta in 2000, 2005, and 2010.

Figure 15. The growth pattern of ISA from 2000 to 2010 in China.

Figure 15. The growth pattern of ISA from 2000 to 2010 in China.

6. Conclusions

Variation in urban areas at large scales has been regarded as a major interface between socio-economic and environmental processes. Gathering accurate and timely information about urban dynamics at national and global scales is especially important for discussing environmental and ecological issues. This paper proposed a new index called NUACI for depicting urban characteristics. It is calculated by integrating DMSP-OLS, MODIS EVI, and NDWI. EVI and NDWI were incorporated into the DMSP-OLS data to extract urban dynamics because of their complementary characteristics among different land covers. DMSP-OLS, VANUI, HSI and NUACI were carried out by means of not only latitudinal transects sampling but also typical cities sampling analysis. The results indicated that NUACI yields the highest performance for reducing pixel saturation and eliminating the blooming effect. Accurate urban information with the highest Kappa coefficient values and best overall accuracy can also be retrieved from the NUACI images by selecting the same sample points across different years.

The NUACI-based regression model was then established to map ISA using the urban fraction derived from TM/ETM+ as reference data. The validation indicated that the NUACI-based regression model has the highest R² value (0.8079), the highest R value (0.8908), and the lowest RMSE value (0.1176). The scatterplots generated with the estimations and the reference data also demonstrated that there is a higher correlation between ISA and NUACI than between ISA and other urban indices. Although there are some variations in the parameters of regression models for different years, the NUACI-based regression models in 2000, 2005, and 2010 also presented higher values for R² and R as well as lower values for RMSE. The ISA results extracted from the NUACI-based regression models indicate that a large proportion of ISA are distributed in the eastern regions, especially in the Pearl River Delta, the Yangtze River Delta, and the Beijing-Tianjin-Tangshan metropolitan areas. The three metropolitan areas show the highest ISA expansion rates during the period from 2000 to 2010.

In summary, the NUACI index proposed in this paper, combining information from multi-source remotely sensed data, can obtain more detailed urban characteristics than single source data can do individually. Our assessments confirm that NUACI has the ability to reduce pixel saturation of NTL, eliminate the blooming effect and provide a more accurate ISA estimation. The NUACI index is easily interpreted and is simple to implement for mapping national urban fractions from only a few samples. Furthermore, the DMSP-OLS and MODIS data used in NUACI can be freely downloaded at the global scale, thus this index is capable of mapping ISA for any large area with limited effort and at low cost in a reasonable amount of time. Furthermore, NPP VIIRS images with newer-version NTL data, with a higher spatial resolution and a wider radiometric detection range will soon become an even more attractive data source for mapping urban areas. We expect that the proposed NUACI will be implemented with NPP VIIRS data to retrieve accurate and timely urban information of even higher resolution at the national and global scales.

We have published the ISA results for China for the years 2000, 2005 and 2010 on the Internet [41], which are available for download and free to use.