Evaluation of Urbanization Dynamics and its Impacts on Surface Heat Islands: A Case Study of Beijing, China

As the capital of China, Beijing has experienced a continued and rapid urbanization process in the past few decades. One of the key environmental impacts of rapid urbanization is the effect of urban heat island (UHI). The objective of this study was to estimate the urbanization indexes of Beijing from 1992 to 2013 based on the stable nighttime light (NTL) data derived from the Defense Meteorological Satellite Program’s Operational Line Scanner System (DMSP/OLS), which has became a widely used remote sensing database after decades of development. The annual average value nighttime light Digital Number (NTL-DN), and the total lit number and urban area proportion within Beijing’s boundary were calculated and compared with social-economic statistics parameters to estimate the correlation between them. Four Landsat thematic mapper (TM) images acquired in 1995 and 2009 were applied to estimate the normalized difference vegetation index (NDVI) and normalized land surface temperature (LSTnor), and spatial correlation analysis was then carried out to investigate the relationship between the urbanization level and NDVI and LSTnor. Our results showed a strong negative linear relationship between the NTL-DN value and NDVI; however, in contrast, a strong positive linear relationship between existed between the NTL-DN value and LSTnor. By conducting a spatial comparison analysis of 1995 and 2009, the vegetation coverage change and surface temperature difference were calculated and compared with the NTL-DN difference. Our result revealed that the regions of fast urbanization resulted in a decrease of NDVI and increase of LSTnor. In addition, choropleth maps showing the spatial pattern of urban heat island zones were produced based on different temperatures, and the analysis result indicated that the spatial distribution of surface temperature was closely related with the NTL-DN and NDVI. These findings are helpful for understanding the urbanization process as well as urban ecology, which both have significant implications for urban planning and minimize the potential environmental impacts of urbanization in Beijing.


Introduction
Currently, more than 50% of the world's population lives in urban areas, with this proportion likely to keep increasing in developing countries. Urbanization is taking place at a spectacular rate worldwide, particularly in China during the past few decades. With continuous rural-urban migration and urban expansion, more megacities are continuously emerging in China [1,2]. Along with rapid urbanization and modernization, some environmental issues, such as water pollution [3], air pollution [4], greenhouse gas emissions [5], and enhanced urban heat islands [6], are arising in remote-sensing data, including the integrated application of DMSP/OLS NTL imagery and Landsat TM products, and estimated the interrelationship between the spatial distribution of urbanization level with the LST nor and NDVI. Furthermore, the comparison analysis was carried out from the DMSP/OLS NTL in 1995 and 2009 to probe into the variation tendency of LST nor and NDVI in different years. In addition, the target region was divided into four categories based on temperature distribution character in each year, and statistical and comparative analysis among the different urban heat island zones were used to explore the correlation between urban heat island effects, and the NTL and NDVI.

Data and Methods
As the capital of China, Beijing is well known worldwide as a famous ancient city and an international metropolis, which is located at 39 • 56' N and 116 • 20' E. Figure 1 shows the geographic location of the study area ( Figure 1a) and administrative districts of the Beijing city ( Figure 1b). As the political, cultural and economic center of the People's Republic of China, Beijing city covers an area of 16,808 km 2 . In terms of topography, two-thirds of Beijing are mountainous areas where the city is surrounded by mountains from all sides, except for the southeast of the city, where a plain slopes slightly to the Bohai Rim. Beijing city is characterized by a warm temperature zone and has a typical continental monsoon climate with four distinct seasons, including a hot and rainy summer and a cold and dry winter.
Remote Sens. 2017, 9,453 3 of 16 urbanization level with the LSTnor and NDVI. Furthermore, the comparison analysis was carried out from the DMSP/OLS NTL in 1995 and 2009 to probe into the variation tendency of LSTnor and NDVI in different years. In addition, the target region was divided into four categories based on temperature distribution character in each year, and statistical and comparative analysis among the different urban heat island zones were used to explore the correlation between urban heat island effects, and the NTL and NDVI.

Data and Methods
As the capital of China, Beijing is well known worldwide as a famous ancient city and an international metropolis, which is located at 39°56' N and 116°20' E. Figure 1 shows the geographic location of the study area ( Figure 1a) and administrative districts of the Beijing city ( Figure 1b). As the political, cultural and economic center of the People's Republic of China, Beijing city covers an area of 16,808 km 2 . In terms of topography, two-thirds of Beijing are mountainous areas where the city is surrounded by mountains from all sides, except for the southeast of the city, where a plain slopes slightly to the Bohai Rim. Beijing city is characterized by a warm temperature zone and has a typical continental monsoon climate with four distinct seasons, including a hot and rainy summer and a cold and dry winter. After 1949 planners were concerned with defining the city as a political and cultural center, the city has been experiencing unprecedented urbanization over the last several decades [27]. The rapid economic growth of Beijing is not only associated with the large-scale migration of population [28], but is also connected with many environmental consequences including UHI effects. Urbanization is a quite complicated phenomenon, with many driving forces responsible for the urbanization of a city. As the capital city of China, the situation in Beijing is more complicated, and it is important to evaluate the transformation processes to support policy making for environmental management and urban planning. Thus, investigating the spatial and temporal urbanization situation and understanding the UHI development situation in Beijing is crucial for the sustainable development of the city in the context of policy making. For the present study, the following methodology was adopted, which involved the collection of social-economic statistics, satellite data collection, preprocessing of the remote sensing imagery, preparation of NDVI maps, retrieval of LSTnor maps, construction of choropleth maps of urban heat island zones and correlation analysis studies. The data sets used in this research are shown in Table 1 below. After 1949 planners were concerned with defining the city as a political and cultural center, the city has been experiencing unprecedented urbanization over the last several decades [27]. The rapid economic growth of Beijing is not only associated with the large-scale migration of population [28], but is also connected with many environmental consequences including UHI effects. Urbanization is a quite complicated phenomenon, with many driving forces responsible for the urbanization of a city. As the capital city of China, the situation in Beijing is more complicated, and it is important to evaluate the transformation processes to support policy making for environmental management and urban planning. Thus, investigating the spatial and temporal urbanization situation and understanding the UHI development situation in Beijing is crucial for the sustainable development of the city in the context of policy making. For the present study, the following methodology was adopted, which involved the collection of social-economic statistics, satellite data collection, preprocessing of the remote sensing imagery, preparation of NDVI maps, retrieval of LST nor maps, construction of choropleth maps of Remote Sens. 2017, 9, 453 4 of 16 urban heat island zones and correlation analysis studies. The data sets used in this research are shown in Table 1 below. The indicators were primarily selected according to the following general selection criterion: (1) choose the most cited indicators; (2) cover the components of urbanization processing and (3) choose the simplest indicators to facilitate indicators. Permanent population and gross domestic product (GDP) were selected as representative indicators of demographic and economic, in addition, these two indicators were the most widely used socio-economic indexes to establish the relationship with the urbanization processes based on DMSP/OLS NTL, which have been estimated by previous researchers [9,10]. Moreover, studies have confirmed that energy consumption and electricity consumption have close relationships with DMSP/OLS NTL based urbanization estimating [13,30]. Besides, popularization and promotion of natural gas in some ways could reflect the development of energy industry structure in the city, there is also little research on establishing relationship with DMSP/OLS NTL, as a result, total consumption of natural gas was selected. Freight traffic plays an important role during the development of a city, length of highways and the total number of civil motor vehicles were chose to reflect the freight traffic development in Beijing. Therefore, the data collection includes seven annual indexes: GDP; permanent population; total energy consumption; total consumption of natural gas; the total electricity consumption; length of highways; and the total number of civil motor vehicles.  [31]. Annual composites were produced for each satellite using the highest-quality of collected data. These images are grid-based compositions with a 0-63 digital number (DN) and a 30 arc-second (approximately 1 km at the equator) spatial resolution for pixels.

Inter-Calibration
Due to the lack of on-board calibration, the NTL-DN among annual DMSP/OLS NTL images were not compatible. To improve the comparability of NLT data in Beijing from the period 1992-2013, a second order regression model proposed by Elvidge et al. [32] was utilized to generate an inter-calibrated DMSP/OLS NTL time series and was matched with the composite of F12 in 1999 to minimize the effects of variation among sensors. On this basis, the NTL series can be adjusted to the same radiometric baseline [32]. As the model was empirically developed using the region with little NTL change as the reference data, by analyzing the socio-economic characteristics based on GDP and built-up area data for Chinese cities from 1992-2013, Yichun City in Heilongjiang Province, China was selected as the reference region. A set of annual composites were compared with F121999 for the city of Yichun using the regression model below. In Equation (1), DN is the original NTL-DN value, DN calibrated is the inter-calibrated, a, b, and c are coefficients.

Intra-Annual Composition and Inter-Annual Correction
To make full use of the data collected from two satellites for the same year and to remove the intra-annual unstable lit pixels, a systematic correction method proposed by Liu et al. [33] was applied as intra-annual composition. The composition produced one intra-annual composite for the year with two NTL images from the same year.
In addition, to make sure the NTL based urban dynamics matched the actual process and to remove the discrepancies in the multi-year dataset, an assumption was proposed that the area would develop continuously, leading to an only grow trend of NTL-DN value over time [31]. The inter-annual correction carried out was applied to each NTL image [33].

Calculation of DMSP/OLS NTL Indicators
Our method for analyzing the dynamics of urbanization processes in Beijing includes three main steps: calculating the annual average NTL-DN; counting the total number of lit pixels; calculation the urban area proportion.

Calculation of the annual average NTL-DN
In this study, we used the annual average NTL-DN as one of three indicators to estimate the annual degree of urbanization in Beijing. The annual average NTL-DN was defined as the annual average NTL-DN value of all pixels located in the administrative districts. By taking the pixels located at the boundaries into consideration, the DN value of each pixel was multiplied by the area that it located. In addition, the calculation was made after the preprocessing of the DMSP/OLS NTL imagery, and can be illustrated by Equation (2). In the equation, DN pixel and Area pixel are the NTL-DN value and the trapezoid area of each pixel, Area sum is the sum of the area within the administrative boundary.
2. Counting the total number of lit pixels The total number of lit pixels explained the overall level of how many pixels were lit within the scope of the study region. The total number of pixels derived from the satellite within Beijing city was 24,991, and after inter-annual correlation, the total number of pixels greater than 0 in each year from each satellite was counted. For the year with two satellite datasets, the pixel was counted sonly if it was lit in either dataset.

Urban area proportion
For the extraction of lit urban area, we adopted the threshold technique developed by Henderson et al. [34]. Furthermore, the optimal thresholds for extraction were determined based on the method provided by Gao et al. [35]. As Beijing is in the Northern Coastal China (NCC) region, the thresholds of NCC were exploited for the extraction.

Preprocessing of Landsat TM Images
Four Landsat 5 TM images were collected on 16 September 1995 (Landsat Scene ID: LT51230321995259HAJ00 and LT51230331995259HAJ00) and 22 September 2009 (Landsat Scene ID: LT51230322009265IKR00 and LT51230332009265IKR00) from the United States Geological Survey (USGS) Earth Explorer website [36]. The Landsat TM images were all acquired under clear sky conditions. The images were further rectified to the Universal Transverse Mercator projection system (datum WGS84, UTM zone N50). The Landsat TM products were composed of seven bands, six of them in the visible and near infrared, with only one band located in the thermal infrared region.
To cover the whole administrative region of Beijing, two images were merged together as a preliminary step. The next step consisted of the calculation of at-sensor spectral radiance in converting image data into a physically meaningful common radiometric scale. Equation (3) was applied to convert the digital numbers for both reflective and thermal bands to at-sensor radiance [37,38]: where L λ is the spectral radiance at the sensor's aperture in W/(m 2 sr µm); Q cal is the quantized calibrated pixel value in DN; Q calmin is the minimum quantized calibrated pixel value corresponding to LMI N λ ; Q calmax is the maximum quantized calibrated pixel value corresponding to LMAX λ ; LMI N λ is the spectral radiance that is scaled to Q calmin in W/(m 2 ·sr·µm); and LMAX λ is the spectral radiance that is scaled to Q calmax in W/(m 2 ·sr·µm).
To reduce the between-scene variability of relatively clear Landsat scenes for all bands, we adopted top-of-atmosphere (TOA) reflectance to correct the cosine effect of solar zenith angles and solar irradiance [37]. The combined surface and atmospheric reflectance was computed as per Equation (4) [38]: where ρ λ is the unitless planetary reflectance; L λ is the spectral radiance at the sensor's aperture in W/(m 2 ·sr·µm); d is the earth-sun distance (in astronomical units); ESUN λ is the mean solar exo-atmospheric irradiance in W/(m 2 ·µm); and θ s is the solar zenith angle in degrees.

Normalized Difference Vegetation Index Calculation
The NDVI index is a measure of the amount and vigor of vegetation at the surface. As vegetation reflects well in the near infrared part of the spectrum, NDVI has become a simple graphical indicator to assess the vegetation coverage of the target. Numerous studies have focused on understanding the relationship between LST and NDVI [26] and concluded that a close relationship between them existed. For this study, to explore the relationship between vegetation coverage and LST and DMSP/OLS based urbanization level, the Landsat image based NDVI was calculated as the reflection of the vegetation coverage. Data supplied by Band 3 and Band 4 can be used to construct this vegetation index according to the following equation:

Estimating Land Surface Temperature
The procedure described by Weng et al. [39] was adopted for the retrieval of LST. The thermal infrared band (Band 6) was utilized to map LST, and the processing was composed of four parts:

1.
First, converting spectral radiance to at-sensor spectral, which is a more physically useful variable, and was calculated by Equation (3); 2. Second, to converting the result of Step 1 to at-sensor brightness temperature, this proceeded under the assumption that the earth's surface is a black body, and using the formula below [38]; where T B is the effective at sensor brightness temperature in Kelvin; L λ is the spectral radiance at the sensor's aperture in W/(m 2 ·sr·µm); K 1 and K 2 are the calibration constant for Landsat TM; K 1 is 607.76 W/(m 2 ·sr·µm); and K 2 is 1260.56 K.

3.
Third, the obtained brightness temperatures were referenced to a black body, which differ from the properties of real objects. Correction for spectral emissivity (ε) is a must and the emissivity corrected surface temperature was computed as follows [40]: where T s is the surface radiant temperature in Kelvin (k); T B is the effective at sensor brightness temperature in Kelvin; λ is the wavelength of emitted radiance, herein, λ = 11.5 µm [41]; Js ; c is the velocity of light 2.998 × 10 8 m/s ; and ε is the surface emissivity which can be calculated from NDVI [42,43], the estimation was shown in Table 2 below.
As the surface radiant temperature is in Kelvin, which differs from the commonly used centigrade, the LST in Celsius was calculated by adding the absolute zero (approximately −273.15 • C) [44].

4.
Finally, to ensure that the images from different years were comparable, normalized land surface temperature was put forward, which quantified the LST range from 0-1. The calculation required the equation below: where LST nor is the normalized land surface temperature; LST min and LST max are the minimum and maximum value obtained within the study region, respectively.

Urbanization Dynamics Process Assessment in Beijing
In this study, the annual average NTL-DN value under Beijing's boundary was calculated using Equation (2). The dynamics of urbanization, which was estimated by DMSP/OLS NTL, are shown in Figure 2a. By preprocessing the images, the calibrated annual average NTL-DN value was calculated and is displayed in red. The average NTL-DN increased 1.56 times, from 11.73 in 1992 to 18.28 in 2013, with an average annual growth rate of 2.16%. Via calibration, we obtained a smooth result of an average NTL-DN trend, which seems more reasonable than the initial NTL data. The total lit number and urban area proportion trends are shown in Figure 2b  3.46 times from 4.86% in 1992 to 16.83% in 2013. Thus, the lit number increased more rapidly in the first decade than the next; however, the urban area proportion showed an opposite trend where it increased faster in the later decade than before.

Correlation Analysis of the NTL Indicators with Social-Economic Variables
The relationship between the NTL indicators and social-economic variables was initially examined through a log-linear regression analysis [45]. Representatively, the Average DN-Ln(GDP), Total Lit Number-Ln(GDP) and Urban Area proportion-Ln(GDP) are shown in Figure 3. The regression analysis determined that the DMSP/OLS based indicators were closely related to the GDP of Beijing during the study period. The correlation coefficient of determination of Ln(GDP) with average NTL-DN value, total lit number and urban area proportion were 0.924, 0.807 and 0.815, respectively, with both showing a close positive correlation. With in-depth analysis through a comparison of the correlation of GDP with NTL indicators, we found that the average NTL-DN had a better correlation than the other two indexes from the approximate increase trends of GDP and Average NTL-DN. Log-linear regression was then applied between the NTL indicators and all collected social-economic parameters from 1992-1993 to estimate their correlation. Table 3 shows the correlation results, and both average NTL-DN and total lit number and urban area proportion were closely correlated with the social-economic variable. The correlation coefficients of determination of annual average NTL-DN value were all above 0.752 with an average value of 0.873. In particular, the correlation of average NTL-DN with the total consumption of natural gas is among the best, with a

Correlation Analysis of the NTL Indicators with Social-Economic Variables
The relationship between the NTL indicators and social-economic variables was initially examined through a log-linear regression analysis [45].

Correlation Analysis of the NTL Indicators with Social-Economic Variables
The relationship between the NTL indicators and social-economic variables was initially examined through a log-linear regression analysis [45]. Representatively, the Average DN-Ln(GDP), Total Lit Number-Ln(GDP) and Urban Area proportion-Ln(GDP) are shown in Figure 3. The regression analysis determined that the DMSP/OLS based indicators were closely related to the GDP of Beijing during the study period. The correlation coefficient of determination of Ln(GDP) with average NTL-DN value, total lit number and urban area proportion were 0.924, 0.807 and 0.815, respectively, with both showing a close positive correlation. With in-depth analysis through a comparison of the correlation of GDP with NTL indicators, we found that the average NTL-DN had a better correlation than the other two indexes from the approximate increase trends of GDP and Average NTL-DN.  Log-linear regression was then applied between the NTL indicators and all collected social-economic parameters from 1992-1993 to estimate their correlation. Table 3 shows the correlation results, and both average NTL-DN and total lit number and urban area proportion were closely correlated with the social-economic variable. The correlation coefficients of determination of annual average NTL-DN value were all above 0.752 with an average value of 0.873. In particular, the correlation of average NTL-DN with the total consumption of natural gas is among the best, with a  Log-linear regression was then applied between the NTL indicators and all collected social-economic parameters from 1992-1993 to estimate their correlation. Table 3 shows the correlation results, and both average NTL-DN and total lit number and urban area proportion were closely correlated with the social-economic variable. The correlation coefficients of determination of annual Remote Sens. 2017, 9, 453 9 of 16 average NTL-DN value were all above 0.752 with an average value of 0.873. In particular, the correlation of average NTL-DN with the total consumption of natural gas is among the best, with a determination coefficient of 0.966. Furthermore, the best correlation of total lit number with social-economics parameters was with the total consumption of natural gas with a determination coefficient of 0.954, while the average correlation determination coefficient was 0.766. For the correlation analysis with urban area proportion, the coefficients of determination were all above 0.674, and held an average value of 0.804, while the maximum of 0.916 was derived from the population. The result illustrates that DMSP/OLS NTL based urbanization indicators are effective techniques to evaluate and monitor development at city level.   Figure 4c,f, respectively). From the spatial distribution in both years, we concluded that the downtown area had a smaller NDVI value and a larger value of LST nor , and that through a visual inspection, a surface urban heat island existed within the city.

Spatial Relationship Analysis of DMSP/OLS NTL and Landsat Images
An analysis based on the DN value was also conducted to explore the spatial correlation between the DMSP/OLS NTL-DN and Landsat based parameters. The calibrated NTL images are shown in  (Figure 5f) were established, and in contrast to the mean NDVI, the LST nor showed a close positive correlation with determination coefficients reaching 0.651 in 1995 and 0.926 in 2009. In both years, the LST nor showed an increasing trend along with the rise of the NTL-DN value, which can be used to draw conclusions on how the urbanization level impacts the temperature distribution. The results indicate that for an area with a higher NTL-DN value (that is more urbanized), the mean LST nor within it would be higher than the surroundings, which provides strong evidence that the surface urban heat island effects are influenced by urban development. and in 2009 (Figure 4e), an obvious horizontal gradient appeared between the urban area and its surroundings in both 1995 and 2009. The value of the NDVI ranges from −0.63-0.79 in 1995, and changes to −0.93-0.75 in 2009, and it can be seen that the NDVI changed significantly due to urban expansion during this period. Furthermore, the LST in 1995 and 2009 were normalized (shown in Figure 4c,f, respectively). From the spatial distribution in both years, we concluded that the downtown area had a smaller NDVI value and a larger value of LSTnor, and that through a visual inspection, a surface urban heat island existed within the city. The results indicate that for an area with a higher NTL-DN value (that is more urbanized), the mean LSTnor within it would be higher than the surroundings, which provides strong evidence that the surface urban heat island effects are influenced by urban development.

Spatial-Temporal Comparison Analysis in 1995 and 2009
To investigate the influence of SUHI driven by urban development, the images of NTL-DN, NDVI and LST nor from two years was composited to obtain the difference between 1995 and 2009. Next, a spatial-temporal comparison analysis was carried out to explore the variation in vegetation coverage situation and the SUHI effects from the urbanization process revealed by NTL images. Figure 6a shows the spatial pattern of the NTL-DN difference between the two years through visual assessment, with the NTL-DN of the previous downtown area remaining unchanged, whereas the surrounding areas had an obvious increase in NTL-DN value, which matched the urban sprawl and expansion in Beijing during this period. With the extraction of the NTL-DN difference, correlation analysis was then carried with using the NDVI difference during the study period (Figure 6b). The NDVI difference showed a negative correlation with the NTL-DN difference, where the determination coefficient reached 0.681, and the correlation was more obvious in the range where theNTL-DN difference had a positive value. The result suggested that the urbanization process would reduce vegetation coverage, given that urban sprawl changed the vegetation area into impervious surfaces. Another correlation analysis was conducted with the NTL-DN difference and LST nor difference (Figure 6c), where the result clearly demonstrated a positive correlation between the two parameters, where the coefficient of determination was 0.664. Areas with a higher value of NTL-DN difference, which represent faster urbanization zones, and the LST nor difference showed a tendency to increase. The faster the development in the area, the higher the increase of LST nor difference. Additionally, the result further illustrated how urbanization affected the SUHI quantitatively over time in the city.

Spatial-Temporal Comparison Analysis in 1995 and 2009
To investigate the influence of SUHI driven by urban development, the images of NTL-DN, NDVI and LSTnor from two years was composited to obtain the difference between 1995 and 2009. Next, a spatial-temporal comparison analysis was carried out to explore the variation in vegetation coverage situation and the SUHI effects from the urbanization process revealed by NTL images. Figure 6a shows the spatial pattern of the NTL-DN difference between the two years through visual assessment, with the NTL-DN of the previous downtown area remaining unchanged, whereas the surrounding areas had an obvious increase in NTL-DN value, which matched the urban sprawl and expansion in Beijing during this period. With the extraction of the NTL-DN difference, correlation analysis was then carried with using the NDVI difference during the study period (Figure 6b). The NDVI difference showed a negative correlation with the NTL-DN difference, where the determination coefficient reached 0.681, and the correlation was more obvious in the range where theNTL-DN difference had a positive value. The result suggested that the urbanization process would reduce vegetation coverage, given that urban sprawl changed the vegetation area into impervious surfaces. Another correlation analysis was conducted with the NTL-DN difference and LSTnor difference (Figure 6c), where the result clearly demonstrated a positive correlation between the two parameters, where the coefficient of determination was 0.664. Areas with a higher value of NTL-DN difference, which represent faster urbanization zones, and the LSTnor difference showed a tendency to increase. The faster the development in the area, the higher the increase of LSTnor difference. Additionally, the result further illustrated how urbanization affected the SUHI quantitatively over time in the city.  Figure 7a,d show the spatial pattern of urban heat island zones in 1995 and 2009, respectively. The average normalized LST was calculated within each NTL pixel to indicate the surface temperature distribution. Next, choropleth maps were produced based on the classification scheme of standard deviation, in which the first-class data values were greater than one standard deviation above the mean, and second-class values were between the mean and one standard deviation above the mean, and so on [39,46]. Thus, the area was divided into four categories: heat island zone,  The average normalized LST was calculated within each NTL pixel to indicate the surface temperature distribution. Next, choropleth maps were produced based on the classification scheme of standard deviation, in which the first-class data values were greater than one standard deviation above the mean, and second-class values were between the mean and one standard deviation above the mean, and so on [39,46]. Thus, the area was divided into four categories: heat island zone, sub-heat island zone, medium temperature zone, and low temperature zone. In both 1995 and 2009, the urban heat island zone was obviously assembling in the city center, and the heat island zone had a large expansion from all directions. Figure 7b,e shows the area proportions and the average LST nor for the four categories, where the percentage of urban heat island zone increased from 9.3% in 1995 to 21.7% in 2009; the sub-heat island zone had no significant change from 29.3% to 28.6%; and both the medium temperature and low temperature zones had a decrease of area proportion, especially the low temperature area which decreased from 44.0%-35.7%. These changes indicated an adverse effect on the thermal environment through the urbanization process from 1995-2009, and that the urban heat island zone occupied a higher proportion than before.

Urban Heat Island Analysis
To better understand the relationship between the NTL-DN, NDVI and UHI zones, the statistics of average NTL-DN and NDVI in different UHI zones were obtained by super-imposing the heat island maps with images of NTL and NDVI. From Figure 7c 34, respectively, from the low temperature zone to the heat island zone. The results illustrated the spatial correlation between the surface temperature and the NTL images. The urban heat island zones located in more developed areas, the sub-heat island and medium temperature zones took second and third place, and the low temperature zones were mainly located in the western and northern parts of the city where most of the area was covered with mountains and forests that have lower NTL-DN values than downtown. Furthermore, in both 1995 and 2009, the low temperature zones had the highest value of NDVI, whereas the urban heat island zones had the lowest. These results showed a positive correlation between the proportion of green space and their cooling effects.
Remote Sens. 2017, 9,453 12 of 16 sub-heat island zone, medium temperature zone, and low temperature zone. In both 1995 and 2009, the urban heat island zone was obviously assembling in the city center, and the heat island zone had a large expansion from all directions. Figure 7b,e shows the area proportions and the average LSTnor for the four categories, where the percentage of urban heat island zone increased from 9.3% in 1995 to 21.7% in 2009; the sub-heat island zone had no significant change from 29.3% to 28.6%; and both the medium temperature and low temperature zones had a decrease of area proportion, especially the low temperature area which decreased from 44.0%-35.7%. These changes indicated an adverse effect on the thermal environment through the urbanization process from 1995-2009, and that the urban heat island zone occupied a higher proportion than before.
To better understand the relationship between the NTL-DN, NDVI and UHI zones, the statistics of average NTL-DN and NDVI in different UHI zones were obtained by super-imposing the heat island maps with images of NTL and NDVI. From Figure 7c The results illustrated the spatial correlation between the surface temperature and the NTL images. The urban heat island zones located in more developed areas, the sub-heat island and medium temperature zones took second and third place, and the low temperature zones were mainly located in the western and northern parts of the city where most of the area was covered with mountains and forests that have lower NTL-DN values than downtown. Furthermore, in both 1995 and 2009, the low temperature zones had the highest value of NDVI, whereas the urban heat island zones had the lowest. These results showed a positive correlation between the proportion of green space and their cooling effects.

Result Analysis
Urban systems do not naturally develop, and are always affected by human activities [2]. As a remote sensing database, which is closely related to anthropogenic activities, the time series DMSP/OLS images provide a consistent and timely measure to estimate social-economic dynamics and changes [7]. In this study, the calibrated annual NTL-DN value from 1992-2013 was used to reflect the average urbanization level of the city in that year, as a global index to evaluate the development level of Beijing. The results show that the human activities observed via satellite were consistent with the official statistical data. The average value of the correlation coefficient of determination between the three remote sensing based indexes with the seven statistical indicators reached 0.873, 0.766 and 0.804, respectively. Compared with the other two indicators, the coefficient of total lit number was a little less, most likely due to the characteristics of urbanization pattern in Beijing. The lit number increased quickly in the few years from 1992 as a fast period of urban sprawl, then slowed down and the urbanization level in the expanded area developed rapidly instead of a continuous outward expansion. In contrast, the social-economic statistics showing a sustained and rapid growth during the period led to less correlation with the total lit number than the other two parameters. The urbanization situation in the study period was confirmed by three fields: the average NTL-DN value (within the administrative boundaries of Beijing) presented an average urbanization situation of the city in each year; the total lit number was an indicator from which the number of pixels lit by human activities could be obtained; and the urban area proportion revealed the cases of urban area sprawl during the study period.
Land surface temperature observations acquired from remote sensing technologies such as MODIS and Landsat were applied to assess the UHI, and to analyze the relationship between surface temperature and land use and land cover (LULC) in urban areas [39]. Nighttime light imagery derived from DMSP/OLS has been widely used for detecting human settlement, dynamic mapping of urban areas and exploring the spatial urban sprawl [8,9]. For the use of these two kinds of databases in combination, Liao et al. [19] found that the DMSP/OLS NTL based energy consumption has a significantly positive correlation with the nighttime SUHI, the results indicated that the anthropogenic heat released from energy consumption is an important contributor to the urban thermal environment [19]. In addition, researchers also found that urban patterns extracted using DMSP/OLS were consistent with those extracted using Landsat products [33]. However, the number of studies estimating the correlation between the Landsat based LST with nighttime imagery from the DMSP/OLS are rare and limited. This paper proposed a new way of integrating different remote sensing data to explore the performance of NDVI and LST by Landsat TM images under different values of NTL-DN from DMSP/OLS imagery. The results clearly showed that the spatial pattern of LST nor and NDVI correlated with the pattern of urbanization level in 1995 and 2009. Moreover, by spatial comparison analysis between the two years, the spatial pattern of LST nor increase and NDVI decrease also proved to be closely related to the change of urbanization level index developed from the nighttime imagery. Furthermore, the urban heat island zone maps presented the spatial pattern of surface temperature distribution, where the results revealed that the urban heat island zone held higher NTL-DN values and lower NDVI than other zones.

Limitations and Potential Improvement
Although the DMSP/OLS time series was successfully adopted to explore urbanization dynamics and analyze the UHI situation in Beijing, there are some sources of uncertainty in this study. First, for the inter-calibration, Yichun City of Heilongjiang Province was selected as the reference region even though the area had little change of NTL during the study period; however, the inappreciable development of the city reduced the accuracy of the calibration. Second, for the calculation of the urban area proportion, we adopted the threshold technique provided by Gao et al. [35], where the thresholds for urban area extraction were determined for several periods rather than every year, which is one of the influence factors for the accuracy of extraction. Third, the calculation of the normalized LST contained limitations. The normalization relied heavily on the extreme values within the region; however, as most of the pixel values ranged in the middle area between the maximum and minimum values, normalization reduced the sensitivity of the analysis. Fourth, in this study, four temporal Landsat time series image data were used to characterize and quantify the spatiotemporal patterns of NDVI and LST nor , and correlation analysis was carried out with the annual value of NTL-DN from DMSP/OLS. As this would have created limitations for the correlation analysis between the annual composites of DMSP/OLS and Landsat images, it was necessary to develop new methods that could be applied effectively to detect the annual value of NDVI and LST nor using time series image data.

Conclusions
This study developed the annual average NTL-DN value, total lit number and urban area proportion to investigate the spatial-temporal urbanization progress in Beijing from 1992-2013. Our results offer quantitative measures for evaluating the urbanization level from three aspects over the period, and provide a new way to explore the urbanization process that is relatively accurate and comprehensive.
From both the DMSP/OLS NTL indicators and the social-economic statistics data, the results showed that Beijing experienced continued and rapid urbanization during the study period. Log-linear regression analysis between the two kinds of data sets showed close correlation with each other, with the average value of determination coefficient of annual average NTL-DN value, total lit number and urban area proportion reached 0.873, 0.766 and 0.804, respectively.
From the distribution map of NDVI and LST nor , we found that the downtown area had less vegetation coverage and a higher LST nor . Furthermore, according to the correlation analysis with the NTL-DN value, our results indicated that both the mean NDVI and mean LST nor were closely correlated with the NTL data. The NDVI and LST nor showed a tendency of decreasing and increasing, respectively, along with the growth of the NTL-DN value. The study proposed a new way of integration using different remote sensing data to explore the urban heat island effects by Landsat TM images and DMSP/OLS imagery. The results provided quantitative measurements to estimate the vegetation coverage situation and surface temperature distribution, which are significant for the eco-environment assessment of the city.
The study also explored the spatial-temporal difference of NDVI and LST nor along with the urbanization process derived from NTL data in 1995 and 2009. Along with the rapid urbanization of the city, the original land use turned into urban land, which is largely composed of an impervious surface [23]. The correlation analysis result between the NTL-DN difference with a change in NDVI and LST nor showed a close relationship, which provided strong evidence of how the urbanization process affects the thermal environment and vegetation coverage over time. In addition, by zonal classification, we also found that the urban heat island zone occupied a higher proportion than before by a growth rate of 12.4%, and that the average NTL-DN and NDVI showed an obvious gradient from the low temperature zone to heat island zone. Our results illustrated the spatial correlation between the surface temperature and the urban development situation, as well as the cooling effects of the green space, which have important implications for urban planning to mitigate the SUHI effects.