An Improved Vegetation Adjusted Nighttime Light Urban Index and Its Application in Quantifying Spatiotemporal Dynamics of Carbon Emissions in China

Nighttime Light (NTL) has been widely used as a proxy of many socio-environmental issues. However, the limited range of sensor radiance of NTL prevents its further application and estimation accuracy. To improve the performance, we developed an improved Vegetation Adjusted Nighttime light Urban Index (VANUI) through fusing multi-year NTL with population density, the Normalized Difference Vegetation Index and water body data and applied it to fine-scaled carbon emission analysis in China. The results proved that our proposed index could reflect more spatial variation of human activities. It is also prominent in reducing the carbon modeling error at the inter-city level and distinguishing the emission heterogeneity at the intra-city level. Between 1995 and 2013, CO2 emissions increased significantly in China, but were distributed unevenly in space with high density emissions mainly located in metropolitan areas and provincial capitals. In addition to Beijing-Tianjin-Hebei, the Yangzi River Delta and the Pearl River Delta, the Shandong Peninsula has become a new emission hotspot that needs special attention in carbon mitigation. The improved VANUI and its application to the carbon emission issue not only broadened our understanding of the spatiotemporal dynamics of fine-scaled CO2 emission, but also provided implications for low-carbon and sustainable development plans.


Introduction
Nighttime Light (NTL), derived from the Defense Meteorological Satellite Program's Operational Linescan System (DMSP/OLS), has become one of the most widely-used data for interpreting some socio-environmental issues that are difficult to understand at a relatively high spatial resolution, such as economic output [1,2], urbanization [3][4][5], poverty [6][7][8], energy consumption and greenhouse gas emissions [9][10][11].Methods and indicators based on NTL are usually labor saving and unbiased compared to those grounded on statistical data.It is especially crucial for increasing our understanding of the spatiotemporal dynamics of the aforementioned issues and thus helping us to make proper and timely decisions toward sustainability [12].
On the one hand, despite the endeavors made by previous studies to improve NTL data, some defects still exist, which may constrain the further application and estimation accuracy.First, because of the limited radiance range of the OLS sensor, pixels with weak radiance (<10 −10 watts/cm 2 /sr/mm), which often exist in developing countries or hamlets, cannot be detected and are often given the value of zero [13].Consequently, the estimations in these unlit pixels might be missing [14,15].Second, the saturation problem in the downtown area and the blooming effect in the suburbs of DMSP/OLS NTL may cause the estimation to deviate from the reality [16][17][18].To deal with these problems, various methods and indices have been developed generally through fusing NTL with other data sources.For example, some used population density to increase the variation of NTL in urban cores and supplemented information of human activities for unlit pixels [14,19].Assuming that vegetation cover and human activities are inversely correlated, some combined Normalized Difference Vegetation Index (NDVI) with NTL to reduce the saturation effect of NTL [20][21][22].Specifically, the Vegetation Adjusted NTL Urban Index (VANUI), which is defined as (1 − NDVI) × NTL, was a simple and effective index for reducing the saturation effect in urban areas [23].It has been employed as a better proxy than the original NTL in studying urban impervious surfaces [24].However, some research also argued that VANUI might result in estimation error in some regions with a diverse natural environment and vegetation cover, such as western cities in China [24,25].Additionally, the investigations using long-term time series VANUI are still not fully addressed.The development of improved methods and indicators that have high accuracy and can effectively address the spatiotemporal dynamics of socio-environmental issues are increasingly needed.
On the other hand, CO 2 emissions and their mitigation have become one of the most crucial challenges for sustainable development faced by all of the countries in the world.Fossil fuel-related CO 2 emissions contributed to more than 70% of the total greenhouse gas emissions, which cause profound impacts on physical, biological and socioeconomic systems, such as global warming, sea level rise, biodiversity loss, agricultural productivity decline and urban heat islands [26][27][28].However, the lack of authentic and reliable statistical data on energy consumption prevents a better investigation of the spatiotemporal dynamics of CO 2 emissions at a finer scale, which is crucial for local governments to adapt to and mitigate climate change efficiently [29][30][31].Taking China, the largest carbon emitter since 2006 [32], as an example, it has set a 60~65% reduction target of CO 2 emissions intensity by 2030 below its 2005 level [33].In order to decompose the mitigation target to local administrative units (for example, cities), it is important to understand the spatial distribution of CO 2 emissions and its change at various scales [34].Many studies have explored the spatiotemporal dynamics at the country or province level based on statistics [35,36].However, the biggest challenge faced by policy makers and researchers is the data availability at the prefectural city or even finer levels since the detailed information on the amount and structures of energy consumption at those levels is not open to the public or fragmented.Luckily, NTL provides a possible solution to fill this gap.Previous studies have proven that it was feasible to estimate CO 2 emissions at the 1-km spatial level using NTL [37].Yet, few studies focused on the potential estimation errors caused by NTL at a fine scale.For example, CO 2 emissions in unlit pixels might be missing; emissions in urban core areas might be underestimated, caused by the saturation problem; while emissions in suburban areas might be overestimated because of the blooming effect.Again, it is quite necessary to develop improved methods and indicators based on NTL to expand our knowledge of carbon emissions.
To this end, we developed an improved VANUI, which combines NTL with time series NDVI, population density and water body data to calibrate the saturation and blooming problems and to reduce the estimation error in the urban core and rural areas and, then, applied it to the carbon emission issue in China, so that the performance of the proposed index can be tested and our understandings of a fine-scaled CO 2 emission dynamics can also be increased.This paper is organized as follows: Section 2 describes the study area, as well as the data we used.Section 3 introduces the methods for developing an improved VANUI, downscaling statistically accounted CO 2 emissions to fine spatial scales, verifying the model accuracy and analyzing the spatiotemporal dynamics of CO 2 emissions.Section 4 analyzes the results.Section 5 discusses the performance of the proposed new index and provides potential uses for further studies, followed by conclusions in Section 6.

Study Area
The jurisdiction units in China are organized under a country-province/municipality-prefectural city-county hierarchy.As shown in Figure 1, detailed energy statistical data for accounting for CO 2 emissions are only available for the whole country, 30 provinces/municipalities (Tibet lacks data) and 30 prefectural cities from 1995.Table 1 provides a summary of the socioeconomic conditions for those 30 cities having detailed energy statistics.However, the same energy inventories for most cities and counties are not available, thus they cannot facilitate a long-term CO 2 emissions study at those scales.

Data Description
Table 2 gives a summary of data used for analysis, which generally include energy statistics and seven kinds of geospatial data.
In order to calculate CO2 emissions from energy consumption in China at provincial and some specific prefectural city level, the amounts of 11 types of energy were collected based on provincial and local energy balance tables from China's Energy Statistical Yearbook [38] and relevant City Statistical Yearbooks of 30 prefecture-level cities as shown in Figure 1.In addition, heat and electricity were also included in the CO2 accounting scopes, which were obtained from the same sources.
The Version 4 DMSP/OLS NTL series record lights on the Earth's surface generated from persistent light in cities, towns and other human settlements and other light from gas flares, fires and illuminated marine vessels.The annual average brightness was presented as six-bit digital numbers (DN) from 0-63 after compositing the highest quality data avoiding interference from sunlit, glare, moonlit, clouds and features from the aurora in a whole year.They spanned −180-180 degrees in longitude and −65-75 degrees in latitude, with 30 arc second grids.
The maximum Moderate-resolution Imaging Spectroradiometer (MODIS) NDVI time series were calculated based on the MOD09GA images from EOS/MODIS using Maximum Value Composites (MVC), after splicing, cutting, projection conversion and unit conversion.The space resolution of the dataset was 500 m × 500 m, and the temporal resolution was a month.Images from January 2000-December 2013 were used in this study.Besides, the Global Inventory Monitoring and Modeling System (GIMMS) NDVI dataset was generated based on the NOAA/AVHRR land dataset after radiometric calibration, geometric correction and cloud elimination, with a spatial resolution of 8 km × 8 km and a temporal resolution of 15 days.Images from January 1995-December 1999 were used as a supplement to fill the data gap between 1995 and 1999.In addition, in order to obtain consistent and comparable NDVI, images from January 2000-December 2006 were employed as samples to capture the relationship between NDVIs from two sensors.In order to match the temporal resolution of MODIS NDVI, monthly GIMMS NDVI data were compiled by the MVC method.
Gridded population density data in China of 1995, 2000, 2005 and 2010 were obtained from the Resources and Environmental Sciences Data Center, Chinese Academy of Sciences (RESDC).The multivariate statistical population spatial model was used for these data based on demographic statistics at a county level and land use/cover data at a spatial resolution of 1 km × 1 km.Besides, urban population density, DEM and the total population were used to ensure that the model results were accurate and reasonable [39].

Data Description
Table 2 gives a summary of data used for analysis, which generally include energy statistics and seven kinds of geospatial data.
In order to calculate CO 2 emissions from energy consumption in China at provincial and some specific prefectural city level, the amounts of 11 types of energy were collected based on provincial and local energy balance tables from China's Energy Statistical Yearbook [38] and relevant City Statistical Yearbooks of 30 prefecture-level cities as shown in Figure 1.In addition, heat and electricity were also included in the CO 2 accounting scopes, which were obtained from the same sources.
The Version 4 DMSP/OLS NTL series record lights on the Earth's surface generated from persistent light in cities, towns and other human settlements and other light from gas flares, fires and illuminated marine vessels.The annual average brightness was presented as six-bit digital numbers (DN) from 0-63 after compositing the highest quality data avoiding interference from sunlit, glare, moonlit, clouds and features from the aurora in a whole year.They spanned −180-180 degrees in longitude and −65-75 degrees in latitude, with 30 arc second grids.
The maximum Moderate-resolution Imaging Spectroradiometer (MODIS) NDVI time series were calculated based on the MOD09GA images from EOS/MODIS using Maximum Value Composites (MVC), after splicing, cutting, projection conversion and unit conversion.The space resolution of the dataset was 500 m × 500 m, and the temporal resolution was a month.Images from January 2000-December 2013 were used in this study.Besides, the Global Inventory Monitoring and Modeling System (GIMMS) NDVI dataset was generated based on the NOAA/AVHRR land dataset after radiometric calibration, geometric correction and cloud elimination, with a spatial resolution of 8 km × 8 km and a temporal resolution of 15 days.Images from January 1995-December 1999 were used as a supplement to fill the data gap between 1995 and 1999.In addition, in order to obtain consistent and comparable NDVI, images from January 2000-December 2006 were employed as samples to capture the relationship between NDVIs from two sensors.In order to match the temporal resolution of MODIS NDVI, monthly GIMMS NDVI data were compiled by the MVC method.
Gridded population density data in China of 1995, 2000, 2005 and 2010 were obtained from the Resources and Environmental Sciences Data Center, Chinese Academy of Sciences (RESDC).The multivariate statistical population spatial model was used for these data based on demographic statistics at a county level and land use/cover data at a spatial resolution of 1 km × 1 km.Besides, urban population density, DEM and the total population were used to ensure that the model results were accurate and reasonable [39].
Moreover, the water body distribution map was used as a mask to eliminate the potential interference caused by low NDVI in water bodies in this paper.All images were projected into the Lambert Azimuthal Equal Area Projection with reference to World Geodetic System 1984 (WGS 84), resampled to a spatial resolution of 1 km × 1 km by the nearest neighbor resampling algorithm and extracted by the multilevel administrative boundaries of China.Besides, four Landsat 8 Operation Land Imager (OLI) images in 2013 covering Beijing, Shanghai, Guangzhou and Urumqi were selected as references to assess the estimation accuracy.

Methods
Figure 2 presents the methodological framework, which includes three major steps.Firstly, an improved VANUI was developed based on the integration of NTL, NDVI, population density and water body data.Secondly, statistically accounted CO 2 emissions at the provincial level were downscaled to a 1-km spatial resolution based on the improved VANUI.Additionally, the accuracy of the prosed new index was validated.Thirdly, the spatiotemporal dynamics of CO 2 emissions at the pixel and prefecture-city level were explored.

An Improved Time-Series VANUI
Because the original time series NTL were obtained from different satellite sensors (F12, F14, F15, F16 and F18), they lack continuity and cannot be used directly for long-term analysis.A systematic calibration approach, which was developed by Liu et al. [40] and included inter-calibration, intra-annual composition and inter-annual series correction processes, was employed in this study to calibrate the original NTL data in China.Besides, the following model proposed by Meng et al. [14] was used to adjust NTL in the lit and unlit areas: where NTL pop is the NTL adjusted by population density; NTLcali is the original night time images after a series of calibration processes proposed by Liu et al. [40]; POP is the population density; 0.34 is the weight for the unlit areas based on the detected electricity access in lit and unlit areas of NTL [41].Specifically, the gridded population density of 1995, 2000, 2005 and 2010 was used to calibrate NTL for the years of 1995~1998, 1999~2003, 2004~2008 and 2009~2013, respectively.We used two datasets of NDVI to satisfy the research period in this paper because MODIS NDVI images are unavailable before 2000.However, the monthly NDVI data from different sensor systems could not be calculated directly because they lack consistency and comparability.Early studies have proven that NDVI from one instrument can be inter-calibrated against another by linear regression models [42].To ensure the consistency of NDVI data from two sensors, we built a linear regression model based on monthly 1-km GIMMS NDVI and MODIS NDVI in China from 2000-2006 (Figure 3) and finally obtained the following model.

An Improved Time-Series VANUI
Because the original time series NTL were obtained from different satellite sensors (F12, F14, F15, F16 and F18), they lack continuity and cannot be used directly for long-term analysis.A systematic calibration approach, which was developed by Liu et al. [40] and included inter-calibration, intra-annual composition and inter-annual series correction processes, was employed in this study to calibrate the original NTL data in China.Besides, the following model proposed by Meng et al. [14] was used to adjust NTL in the lit and unlit areas: where NTL pop is the NTL adjusted by population density; NTL cali is the original night time images after a series of calibration processes proposed by Liu et al. [40]; POP is the population density; 0.34 is the weight for the unlit areas based on the detected electricity access in lit and unlit areas of NTL [41].Specifically, the gridded population density of 1995, 2000, 2005 and 2010 was used to calibrate NTL for the years of 1995~1998, 1999~2003, 2004~2008 and 2009~2013, respectively.We used two datasets of NDVI to satisfy the research period in this paper because MODIS NDVI images are unavailable before 2000.However, the monthly NDVI data from different sensor systems could not be calculated directly because they lack consistency and comparability.Early studies have proven that NDVI from one instrument can be inter-calibrated against another by linear regression models [42].To ensure the consistency of NDVI data from two sensors, we built a linear regression model based on monthly 1-km GIMMS NDVI and MODIS NDVI in China from 2000-2006 (Figure 3) and finally obtained the following model.
where NDVI MODIS and NDVI GIMMS represent monthly 1-km MODIS NDVI and GIMMS NDVI datasets, respectively.The annual mean NDVI time series from 1995-2013 were then calculated.However, unlike MVC, which could reduce the impacts of cloud and shadow of NDVI, the mean NDVI can reduce seasonal sensitivity and fluctuation and, thus, is more stable than the max NDVI [23].Negative NDVI values were removed because they indicate deserts, glaciers and other areas with few human activities [43,44].
The improved VANUI in China from 1995-2013 was then computed using the following model [23]: where NDVI is the annual mean NDVI; NTL nor is the normalized value of the NTL adjusted by population density ranging from 0-1, which can be calculated by: where NTL max and NTL min are the maximum and minimum values of NTL pop from 1995-2013.Similar to Liu's approach to improve the annual consistence of DMSP/OLS NTL during a long time period [40], the annual improved VANUI should be calibrated by Equation ( 5) on the basis of the following assumptions: (1) the pixel level adjusted VANUI grew continuously in China during the last 19 years, and the value in a later year should be equal to or greater than that in the previous year; (2) the adjusted VANUI in 2013 was the largest and closest to reality in the time series.We chose the adjusted VANUI in 2013 instead of VANUI in 1995 as the reference data because VANUI in 2013 was calculated based on MODIS NDVI, which have higher quality than those in 1995.The annual mean NDVI time series from 1995-2013 were then calculated.However, unlike MVC, which could reduce the impacts of cloud and shadow of NDVI, the mean NDVI can reduce seasonal sensitivity and fluctuation and, thus, is more stable than the max NDVI [23].Negative NDVI values were removed because they indicate deserts, glaciers and other areas with few human activities [43,44].
The improved VANUI in China from 1995-2013 was then computed using the following model [23]: where NDVI is the annual mean NDVI; NTL nor is the normalized value of the NTL adjusted by population density ranging from 0-1, which can be calculated by: where NTL max and NTL min are the maximum and minimum values of NTL pop from 1995-2013.Similar to Liu's approach to improve the annual consistence of DMSP/OLS NTL during a long time period [40], the annual improved VANUI should be calibrated by Equation ( 5) on the basis of the following assumptions: (1) the pixel level adjusted VANUI grew continuously in China during the last 19 years, and the value in a later year should be equal to or greater than that in the previous year; (2) the adjusted VANUI in 2013 was the largest and closest to reality in the time series.We chose the adjusted VANUI in 2013 instead of VANUI in 1995 as the reference data because VANUI in 2013 was calculated based on MODIS NDVI, which have higher quality than those in 1995.( where t = 1, 2, . . ., 18; T is the latest year of the study period (that is 2013 in this research); VANUI T−t cali is the VANUI that needs to be calibrated in the year of T − t (i.e., from 2012-1995); VANUI T−t adj is the improved VANUI from 2012-1995 calculated before; VANUI T−t+1 adj is the improved VANUI in the next year of VANUI T−t adj .

Estimating CO 2 Emissions Using Improved VANUI
First, CO 2 emissions from 30 provinces are accounted so as to provide training samples for further downscaling analysis.The guideline proposed by the Intergovernmental Panel on Climate Change (IPCC) was adopted as a standard method for carbon accounting [45]: where SCE pro represents the statistical CO 2 emissions of the sample province; the subscript k denotes the energy type; FC is total fuel consumption; EF is the effective CO 2 emission factor; LC is the low calorific value of each fuel; CC represents carbon content; CR is the carbon oxidation ratio; 44/12 is the molecular weight ratio of carbon dioxide to carbon.The EF of electricity was the mean annual EF of seven power sub-grids in China from 2007-2013, provided by National Development and Reform Commission (NDRC) [46].Others were compiled from the Provincial Greenhouse Gas List published by NDRC [47].Table 3 provides the parameters for emissions accounting.Second, a linear regression model, as shown in Equation (7), was employed to capture the relationship between statistically accounted CO 2 emissions from sample provinces and improved VANUI: where VANUI pro denotes the total improved VANUI value in corresponding provinces; α and β are the intercept and coefficient, which need to be estimated.In this study, 19 regression models from 1995-2013 were built to reduce the temporal bias in estimating CO 2 emissions.As shown in Table 4, all parameters were significant at the 5% or 1% level based on the t-test, with reasonable values of the coefficient of determination R 2 (R 2 > 0.75), suggesting that the improved VANUI is an appropriate proxy of CO 2 emissions.One possible reason for the increasing R 2 after 2000 might be that the quality of VANUI calculated by MODIS NDVI after 2000 is higher than those calculated by GIMMS NDVI because of the higher spatial resolution of the original MODIS NDVI (500 m).Assuming the relationship between CO 2 emissions and improved VANUI at the pixel level was consistent with that at the sample province level, the following model was used to downscale provincial CO 2 emissions to the pixel level: where ĈE pix represents the CO 2 emissions estimated by the improved VANUI at the pixel level; N pro is the total lit pixels in the province; VANUI pix is the value of improved VANUI in the pixel; others are the same as for Equation (7).Furthermore, to limit the estimation residuals within a province, the estimated pixel-level CO 2 emissions by Equation (8) were corrected by the following equation: where CE pix is the modeled pixel-level CO 2 emissions after correction.ĈE pro is the sum of pixel-level estimated CO 2 emissions ( ĈE pix ) within each provincial boundary.It was noted that since Tibet Province lacks energy statistical data, pixel-level CO 2 emissions within Tibet could not be corrected through Equation ( 9), but directly use the estimations by Equation (8).

Assessing the Accuracy of Improved VANUI in Modeling CO 2 Emissions
Due to the lack of high-resolution CO 2 emissions for a reference, we used two assessments to evaluate the performance of improved VANUI in modeling CO 2 emissions.First, the root mean square error (RMSE) and mean relative error (MRE) were used to describe the difference between the modeled CO 2 emissions and statistically accounted emissions from 30 prefecture cities that have detailed energy statistics (the accounting method was the same as that at the provincial level).Second, since CO 2 emissions from impervious areas are usually higher than those from dense vegetation covered areas [36], Landsat 8 OLI-TIRS images of four selected cities in China, including Beijing and Shanghai (typical mega cities that usually have the NTL saturation phenomenon), Guangzhou and Urumqi (typical cities with distinct natural environment and vegetation covers) were used as references to assess whether the improved VANUI could reflect this emission difference within the same city.

Analyzing Spatiotemporal Patterns of Prefecture-City Level CO 2 Emissions
Indicators including Global Moran's I, and Anselin Local Moran's I (Local indicators of spatial association, LISA) were used to explore the spatiotemporal patterns of prefecture-city level CO 2 emissions.Specifically, Global Moran's I measures the spatial autocorrelation of prefectural CO 2 emissions nationwide [48].
where n is the total number of prefectural cities; i and j are two different cities; w ij is the spatial weights matrix, which was decided by the criteria of the first-order queen contiguity in this study; CE i and CE j are the modeled CO 2 emissions for city i and j, respectively; CE is the average emissions of the whole study area.With the value ranging from −1-1, Global Moran's I evaluates whether the spatial distribution of CO 2 emissions is random (=0), clustered (>0) or dispersed (<0).LISA shows the spatial dependence and heterogeneity among prefectural cities [33].
where Z i is the deviation of modeled emissions for city i from the average and Z o i is the spatial lag for city i, which reflects the weighted average of CO 2 emissions from the neighboring city j, and can be calculated by n ∑ j=1 w ij Z j .
Based on whether both Z i and Z o i are higher than 0 or not, cities could be classified into four cluster types.They are high-high clusters (Z i > 0 and Z o i > 0), high-low clusters (Z i > 0 and Z o i < 0), low-high clusters (Z i < 0 and Z o i > 0) and low-low clusters (Z i < 0 and Z o i < 0).The LISA results can be used to provide visual information of local instability in spatial autocorrelation [49].

Performance Test of the Improved VANUI
To evaluate the capability of the improved VANUI in tackling the saturation and blooming problems, we compared the normalized value of original NTL, population-adjusted NTL (calculated by Equations ( 1) and ( 4)), original VANUI and improved VANUI in a line transect in the Beijing-Tianjin-Tangshan metropolitan area.It was found that the values of the original NTL were obviously saturated in the urban center of Beijing (Figure 4b).Although population-adjusted NTL could help to reduce the saturation, the improvements in some urbanized areas were still not as much as expected.For instance, the pixel range of population-adjusted NTL in Figure 4c was still close to one.Notably, both the original and improved VANUI captured more spatial details in urban center areas with their values much lower than other two indices.In suburban areas, the capability of original NTL and original VANUI in describing the human activity variation was weak.For example, the values in the 0-20 pixel range in Figure 4b,d were close to zero and had little change.Large fluctuations in population-adjusted NTL and improved VANUI values could be detected (for example, the pixel range of 12-13 in Figure 4c,e), which suggested that they could supplement more information of the human activity differences.In sum, the improved VANUI can reflect human activities more accurately, which is crucial for interpreting socio-environmental issues.

Accuracy Assessment of Improved VANUI in Modeling CO2 Emissions
Using the same methods described in Section 3.2, CO2 emissions were additionally modeled based on the original NTL, the population-adjusted NTL and the original VANUI and compared to that based on the improved VANUI so that the accuracy of improved VANUI in modeling CO2 emissions can be assessed.The comparisons were conducted at the following two levels.
First, at the inter-city level, as shown in Figure 5, statistically accounted CO2 emissions from 30 prefecture cities that have detailed energy statistics between 1995 and 2013 (in total, 170 samples) were accounted as references by the same method described in Section 3.2.It was found that the R 2 in the correlation between statistical CO2 emissions and modeled emissions derived from the improved VANUI was the same as that in the original VANUI's case (R 2 = 0.64), but was higher than

Accuracy Assessment of Improved VANUI in Modeling CO 2 Emissions
Using the same methods described in Section 3.2, CO 2 emissions were additionally modeled based on the original NTL, the population-adjusted NTL and the original VANUI and compared to that based on the improved VANUI so that the accuracy of improved VANUI in modeling CO 2 emissions can be assessed.The comparisons were conducted at the following two levels.
First, at the inter-city level, as shown in Figure 5, statistically accounted CO 2 emissions from 30 prefecture cities that have detailed energy statistics between 1995 and 2013 (in total, 170 samples) were accounted as references by the same method described in Section 3.2.It was found that the R 2 in the correlation between statistical CO 2 emissions and modeled emissions derived from the improved VANUI was the same as that in the original VANUI's case (R 2 = 0.64), but was higher than that in the original NTL's case (R 2 = 0.59) and the population-adjusted NTL's case (R 2 = 0.52).More importantly, the modeled CO 2 emissions using the improved VANUI owned the lowest RMSE and MAE compared to the others, which suggested that improved VANUI was more accurate in modeling CO 2 emissions at the inter-city level.
that in the original NTL's case (R 2 = 0.59) and the population-adjusted NTL's case (R 2 = 0.52).More importantly, the modeled CO2 emissions using the improved VANUI owned the lowest RMSE and MAE compared to the others, which suggested that improved VANUI was more accurate in modeling CO2 emissions at the inter-city level.Second, at the intra-city level, Figure 6 shows the modeled CO2 emissions (CE) in Beijing, Shanghai, Guangzhou and Urumqi based on four indices.Landsat 8 OLI images in 2013 were used for visual references.The modeled results varied significantly in the four cases.Obviously, the variation of emission density using improved VANUI was the largest with the maximum value around 0.30 million tons/km 2 , which was 4.2-, 1.8-and 2.2-times larger than that in the original NTL, population-adjusted NTL and original VANUI, respectively.This result implied that improved VANUI could strengthen the heterogeneity at the intra-city level.In addition, it was found that CO2 emissions modeled by improved VANUI could reduce the saturation effect and estimation deviation significantly.Due to the blooming problem, emissions based on the original NTL in urban parks and suburbs were usually overestimated, while being underestimated in sub-centers.In the population-adjusted NTL's case, emissions in suburban regions were much lower than those in urban cores, but the saturation phenomenon at the center area of mega cities was still serious.As for the original VANUI-based emissions, though VANUI can capture more details in urban cores, its performance in suburban areas with less vegetation cover (for example, suburban Urumqi) was not satisfactory.Moreover, emissions in urbanized areas with low NDVI, but high NTL might be overestimated (for example, river zones in Shanghai and Guangzhou).By contrast, the relatively low emissions in urban parks and high emissions in urban cores and sub-centers have been successfully distinguished through improved VANUI.Second, at the intra-city level, Figure 6 shows the modeled CO 2 emissions (CE) in Beijing, Shanghai, Guangzhou and Urumqi based on four indices.Landsat 8 OLI images in 2013 were used for visual references.The modeled results varied significantly in the four cases.Obviously, the variation of emission density using improved VANUI was the largest with the maximum value around 0.30 million tons/km 2 , which was 4.2-, 1.8-and 2.2-times larger than that in the original NTL, population-adjusted NTL and original VANUI, respectively.This result implied that improved VANUI could strengthen the heterogeneity at the intra-city level.In addition, it was found that CO 2 emissions modeled by improved VANUI could reduce the saturation effect and estimation deviation significantly.Due to the blooming problem, emissions based on the original NTL in urban parks and suburbs were usually overestimated, while being underestimated in sub-centers.In the population-adjusted NTL's case, emissions in suburban regions were much lower than those in urban cores, but the saturation phenomenon at the center area of mega cities was still serious.As for the original VANUI-based emissions, though VANUI can capture more details in urban cores, its performance in suburban areas with less vegetation cover (for example, suburban Urumqi) was not satisfactory.Moreover, emissions in urbanized areas with low NDVI, but high NTL might be overestimated (for example, river zones in Shanghai and Guangzhou).By contrast, the relatively low emissions in urban parks and high emissions in urban cores and sub-centers have been successfully distinguished through improved VANUI.

Spatiotemporal Dynamics CO2 Emissions in China
Figure 7 shows the spatial pattern of the modeled pixel-level CO2 emissions in China from 1995-2013.It is found that CO2 emissions mainly concentrated in Eastern and Southeastern regions.Visually, the CO2 emission density in metropolitan areas and provincial capital cities was much higher than that in small cities and rural areas.The highly urbanized areas including Beijing-Tianjin-Hebei (BTH), the Yangtze River Delta (YRD) and the Pearl River Delta (PRD) witnessed a sharp increase in emissions, where the emission density has more than doubled from 0.14-0.30 million tons/km 2 between 1995 and 2013.
By aggregating the pixel level CO2 emissions to prefecture cities and analyzing the spatiotemporal characteristics, it was noted that the results of Global Moran's I were significantly positive during the study period, which indicated the existence of spatial agglomeration of CO2 emissions in Chinese cities (Figure 8).In other words, cities with similar emissions tended to cluster in space.However, the degree of agglomeration underwent a decrease from 0.31 in 1995 to 0.20 in 1999, then an increase to 0.36 in 2007 and, finally, became relatively stable around 0.35 from 2008-2013.The differences in regional development strategy and socioeconomic conditions might be one of the main reasons for the changing level of agglomeration.For example, the quick catching up of the economy in inland regions of China due to a series of favorable development strategies and plans, such as the "Grand Western Development Program" since 2000 and the "Rise of Central China Plan" since 2004, may explain the growth of the agglomeration degree in the beginning of the 2000s.Moreover, great efforts were made by the Chinese government to encourage energy

Spatiotemporal Dynamics CO 2 Emissions in China
Figure 7 shows the spatial pattern of the modeled pixel-level CO 2 emissions in China from 1995-2013.It is found that CO 2 emissions mainly concentrated in Eastern and Southeastern regions.Visually, the CO 2 emission density in metropolitan areas and provincial capital cities was much higher than that in small cities and rural areas.The highly urbanized areas including Beijing-Tianjin-Hebei (BTH), the Yangtze River Delta (YRD) and the Pearl River Delta (PRD) witnessed a sharp increase in emissions, where the emission density has more than doubled from 0.14-0.30 million tons/km 2 between 1995 and 2013.
By aggregating the pixel level CO 2 emissions to prefecture cities and analyzing the spatiotemporal characteristics, it was noted that the results of Global Moran's I were significantly positive during the study period, which indicated the existence of spatial agglomeration of CO 2 emissions in Chinese cities (Figure 8).In other words, cities with similar emissions tended to cluster in space.However, the degree of agglomeration underwent a decrease from 0.31 in 1995 to 0.20 in 1999, then an increase to 0.36 in 2007 and, finally, became relatively stable around 0.35 from 2008-2013.The differences in regional development strategy and socioeconomic conditions might be one of the main reasons for the changing level of agglomeration.For example, the quick catching up of the economy in inland regions of China due to a series of favorable development strategies and plans, such as the "Grand Western Development Program" since 2000 and the "Rise of Central China Plan" since 2004, may explain the growth of the agglomeration degree in the beginning of the 2000s.Moreover, great efforts were made by the Chinese government to encourage energy conservation and emissions reduction at the country level, such as the implementation of the "Medium and Long-Term Energy Conservation Plan" launched in 2004, the "11th Five Year Plan (2006-2010)" and the "12th Five Year Plan (2011-2015)", which could be the important reasons why agglomeration of CO 2 emissions kept stable after 2007 [50,51].Figure 9 shows the results of LISA for 343 prefecture cities.Based on four types of clusters, the hot spots of carbon emissions and spatially-heterogeneous areas can be identified.Most of the clusters are spatially coherent.BTH, YRD and PRD belonged to the high-high clusters, suggesting that these areas were hot emission spots that need special attention in carbon mitigation.Notably, because of the rapid development of the Bohai Economic Rim in recent years, Shandong Peninsula has replaced the middle south Liaoning, a traditional heavy industrial region in Northeastern China, and became a new cluster with large emissions.The low-low clusters mainly located in the western region and have become less extensive.Most of these clusters were located in Tibet, Qinghai and Gansu provinces, which have less population and carbon-intensive industries.Chongqing and Nanning were labeled as high-low clusters with more CO2 emissions than their neighboring cities.On the contrary, Chengde, a city next to Beijing, emitted much less CO2 than its neighbors did.

Comparison with Previous Studies
By comparing our modeled CO 2 emissions in some selected cities and years with those from previous studies, most of which are based on statistical data and methods, the robustness and reliability of our proposed methods can be validated.As shown in Figure 10, there were big differences in each study even for the same city based on the same accounting framework.Taking Tianjin as an example, our estimate results were smaller than those from Ref-1, but larger than those of Ref-5.The different accounting scopes and methods might be the main reason for these differences, as Ref-1 accounted the carbon emissions from industrial processes while we did not, and we considered the emissions from electricity and heating, while Ref-5 ignored them.Regardless, the changing trend and scale of all of the results are similar to each other, which suggested the validity of our estimation based on the improved VANUI.

Comparison with Previous Studies
By comparing our modeled CO2 emissions in some selected cities and years with those from previous studies, most of which are based on statistical data and methods, the robustness and reliability of our proposed methods can be validated.As shown in Figure 10, there were big differences in each study even for the same city based on the same accounting framework.Taking Tianjin as an example, our estimate results were smaller than those from Ref-1, but larger than those of Ref-5.The different accounting scopes and methods might be the main reason for these differences, as Ref-1 accounted the carbon emissions from industrial processes while we did not, and we considered the emissions from electricity and heating, while Ref-5 ignored them.Regardless, the changing trend and scale of all of the results are similar to each other, which suggested the validity of our estimation based on the improved VANUI.

Limitations and Potential Uses
Though the proposed methods are proven to be effective and accurate in analyzing the carbon emission issue, there remains several limitations that could be further improved in the future study.First, we did not take into account the uncertainties from emission factors and energy consumption data when accounting CO2 emissions based on statistics, which could vary from ±5-±10% [57,58].It was also pointed out that the poor quality of energy statistics would cause significant discrepancy between national and provincial emissions in China [59].All of these will influence the model accuracy since the statistically accounted CO2 emissions are among the fundamental bases of our methods.
Second, although the estimation error of improved VANUI was the smallest, there still exist big gaps between modeled and statistically accounted emissions for some cities (as shown in Figure 5).

Limitations and Potential Uses
Though the proposed methods are proven to be effective and accurate in analyzing the carbon emission issue, there remains several limitations that could be further improved in the future study.First, we did not take into account the uncertainties from emission factors and energy consumption data when accounting CO 2 emissions based on statistics, which could vary from ±5-±10% [57,58].It was also pointed out that the poor quality of energy statistics would cause significant discrepancy between national and provincial emissions in China [59].All of these will influence the model accuracy since the statistically accounted CO 2 emissions are among the fundamental bases of our methods.
Second, although the estimation error of improved VANUI was the smallest, there still exist big gaps between modeled and statistically accounted emissions for some cities (as shown in Figure 5).
More efforts in improving the data quality and integrating more variables in model algorithms should be made.Specifically, since CO 2 emissions are not only driven by demographic factors, economic level and structure, local biophysical conditions such as temperature, elevation and urban form are also reported as important determinants [60][61][62].In the future, with the emergence and development of big data science, it might be possible to obtain detailed spatial data of the aforementioned variables through the Internet, social media, commercial data companies, and so on, for enhancing the performance of NTL.
Owing to the global coverage of NTL and NDVI and the increasing availability of spatial data for population and land use and cover in the world (for example, Landscan's population data with 1-km resolution [63] and Global Land Project's land use and cover maps with 30-m resolution [64], our proposed methods could also be applied to other countries and regions and to other socio-environmental issues, such as steel stocks' distribution in human society and residential energy consumption.

Conclusions
By improving the original VANUI via the integration of NTL with NDVI, population density and water body distribution, this paper proposed a new index to enhance the performance of DMSP/OLS NTL, so that it can be better applied in analyzing socio-environmental issues.China was selected as a case to validate the accuracy of improved VANUI in modeling fine-scaled CO 2 emissions and to increase our understandings of its spatiotemporal dynamics from 1995-2013.
The findings can be concluded as follows.First, the improved VANUI can supplement more information of human activities in areas with weak radiance, which cannot be detected by the original NTL.Second, it can reduce estimation errors in downtown areas and suburbs caused by the saturation and blooming problems of NTL.Third, the improved VANUI shows great potential to be further applied as a proxy of many socio-environmental and socioeconomic issues, especially in regions that are lack of sufficient and reliable statistical data.Fourth, with the aid of improved VANUI, it is found that carbon emission density in metropolitan areas and provincial capitals was much higher than in small cities and rural areas.Global Moran's I confirmed the existence of the spatially-agglomerated distribution of CO 2 emissions in China.In addition to the BTH, YRD and PRD regions, Shandong Peninsula replaced the middle south Liaoning and became a new emission hotspot.These findings may increase our understanding of the spatiotemporal dynamics of CO 2 emissions and provide decision-making support to design sustainable development plans and low-carbon policies.

Figure 1 .
Figure 1.Study area and cities with detailed energy statistics.

Figure 1 .
Figure 1.Study area and cities with detailed energy statistics.

Figure 2 .
Figure 2. Methodological framework of this study.

Figure 2 .
Figure 2. Methodological framework of this study.

Figure 4 .
Figure 4. Performance test of improved VANUI.(a) A false color composite of normalized NTL, NDVI and improved VANUI in red, green and blue; normalized values of (b) NTL, (c) population-adjusted NTL, (d) original VANUI and (e) improved VANUI in a line transect in the Beijing-Tianjin-Tangshan metropolitan area in 2013.

Figure 4 .
Figure 4. Performance test of improved VANUI.(a) A false color composite of normalized NTL, NDVI and improved VANUI in red, green and blue; normalized values of (b) NTL, (c) population-adjusted NTL, (d) original VANUI and (e) improved VANUI in a line transect in the Beijing-Tianjin-Tangshan metropolitan area in 2013.

Figure 5 .
Figure 5. Correlation analysis between modeled CO2 emissions and statistical CO2 emissions for 30 prefecture cities in 1995-2013.(a) Original NTL; (b) population-adjusted NTL; (c) original VANUI; (d) improved VANUI.RMSE and MRE are the root mean square error and mean relative error, respectively.Mt = million tons.

Figure 5 .
Figure 5. Correlation analysis between modeled CO 2 emissions and statistical CO 2 emissions for 30 prefecture cities in 1995-2013.(a) Original NTL; (b) population-adjusted NTL; (c) original VANUI; (d) improved VANUI.RMSE and MRE are the root mean square error and mean relative error, respectively.Mt = million tons.

Figure 6 .
Figure 6.Comparison of modeled CO2 emissions (CE) based on four indices in urban areas of selected cities with the Landsat 8 images in 2013 as references.

Figure 6 .
Figure 6.Comparison of modeled CO 2 emissions (CE) based on four indices in urban areas of selected cities with the Landsat 8 images in 2013 as references.
Remote Sens. 2017, 9, 829 14 of 20 conservation and emissions reduction at the country level, such as the implementation of the "Medium and Long-Term Energy Conservation Plan" launched in 2004, the "11th Five Year Plan (2006-2010)" and the "12th Five Year Plan (2011-2015)", which could be the important reasons why agglomeration of CO2 emissions kept stable after 2007 [50,51].

Figure 9
Figure9shows the results of LISA for 343 prefecture cities.Based on four types of clusters, the hot spots of carbon emissions and spatially-heterogeneous areas can be identified.Most of the clusters are spatially coherent.BTH, YRD and PRD belonged to the high-high clusters, suggesting that these areas were hot emission spots that need special attention in carbon mitigation.Notably, because of the rapid development of the Bohai Economic Rim in recent years, Shandong Peninsula has replaced the middle south Liaoning, a traditional heavy industrial region in Northeastern China, and became a new cluster with large emissions.The low-low clusters mainly located in the western region and have become less extensive.Most of these clusters were located in Tibet, Qinghai and Gansu provinces, which have less population and carbon-intensive industries.Chongqing and Nanning were labeled as high-low clusters with more CO2 emissions than their neighboring cities.On the contrary, Chengde, a city next to Beijing, emitted much less CO2 than its neighbors did.

Figure 9 .
Figure 9. Spatial dependence and heterogeneity of prefecture city level CO2 emissions.(a-d) are the spatial clusters with a 95% confidence level based on Local Moran's I.

Figure 9 20 Figure 8 .
Figure9shows the results of LISA for 343 prefecture cities.Based on four types of clusters, the hot spots of carbon emissions and spatially-heterogeneous areas can be identified.Most of the clusters are spatially coherent.BTH, YRD and PRD belonged to the high-high clusters, suggesting that these areas were hot emission spots that need special attention in carbon mitigation.Notably, because of the rapid development of the Bohai Economic Rim in recent years, Shandong Peninsula has replaced the middle south Liaoning, a traditional heavy industrial region in Northeastern China, and became a new cluster with large emissions.The low-low clusters mainly located in the western region and have become less extensive.Most of these clusters were located in Tibet, Qinghai and Gansu provinces, which have less population and carbon-intensive industries.Chongqing and Nanning were labeled as high-low clusters with more CO 2 emissions than their neighboring cities.On the contrary, Chengde, a city next to Beijing, emitted much less CO 2 than its neighbors did.

Figure 9 .
Figure 9. Spatial dependence and heterogeneity of prefecture city level CO2 emissions.(a-d) are the spatial clusters with a 95% confidence level based on Local Moran's I.

Figure 9 .
Figure 9. Spatial dependence and heterogeneity of prefecture city level CO 2 emissions.(a-d) are the spatial clusters with a 95% confidence level based on Local Moran's I.

Table 1 .
A summary of administrative area, population, population density and gross domestic product (GDP) for the prefectural cities with detailed energy statistics in 2013.

Table 2 .
Description of data used in this study.

Table 3 .
Parameters for accounting CO 2 emissions from energy consumption in China.CR, carbon oxidation ratio.

Table 4 .
Results of statistical CO 2 emissions and improved VANUI linear regressions.