Modeling Spatiotemporal Population Changes by Integrating DMSP-OLS and NPP-VIIRS Nighttime Light Data in Chongqing, China

The sustained growth of non-farm wages has led to large-scale migration of rural population to cities in China, especially in mountainous areas. It is of great significance to study the spatial and temporal pattern of population migration mentioned above for guiding population spatial optimization and the effective supply of public services in the mountainous areas. Here, we determined the spatiotemporal evolution of population in the Chongqing municipality of China from 2000–2018 by employing multi-period spatial distribution data, including nighttime light (NTL) data from the Defense Meteorological Satellite Program’s Operational Linescan System (DMSP-OLS) and the Suomi National Polar-orbiting Partnership Visible Infrared Imaging Radiometer Suite (NPP-VIIRS). There was a power function relationship between the two datasets at the pixel scale, with a mean relative error of NTL integration of 8.19%, 4.78% less than achieved by a previous study at the provincial scale. The spatial simulations of population distribution achieved a mean relative error of 26.98%, improved the simulation accuracy for mountainous population by nearly 20% and confirmed the feasibility of this method in Chongqing. During the study period, the spatial distribution of Chongqing’s population has increased in the west and decreased in the east, while also increased in low-altitude areas and decreased in medium-high altitude areas. Population agglomeration was common in all of districts and counties and the population density of central urban areas and its surrounding areas significantly increased, while that of non-urban areas such as northeast Chongqing significantly decreased.


Introduction
Urban-rural migration is a major issue affecting the sustainable development of society, while the spatial distribution of population is a core focus of research in population geography [1]. Driven by economic globalization, developing countries occupy an increasing share of the world economy and the world's economic center continues to move to Asia [2][3][4]. As the largest developing country, China has experienced an unprecedented growth rate over the past 30 years. The urbanization rate has increased from 26% to 58% and the growth rate is about 2.7 times the world average (World Bank). China's rapidly developing social economy and ongoing urbanization has resulted in the relocation and reorganization of urban and rural populations [5][6][7][8] as reflected in the continuous growth of the former and substantial reductions in the latter's labor force. According to the National Bureau of Statistics of China (NBSC), the country's urban population has increased by an average of 21 million per year since 2000. In contrast, the agricultural labor force has decreased by 11 million per year [9] and the rural population has decreased by~30.2%, from 808 million in 2000 to 564 million in 2018 (NBSC). It is worth noting that population migration from mountainous areas has been particularly significant [9]. The process of urban-rural migration results in the redistribution of production factors such as capital, which will impact on the ecosystem and social economy, with contradiction between resources, the environment and population changing accordingly [10][11][12]. The rural population structure has also changed (including age, gender and number), which has changed the land use pattern and human activities radius, thereby affecting the construction and restoration of rural ecological civilization [6, [13][14][15]. Therefore, mapping and estimating the spatial distribution of populations can provide scientific support for developing regionally sustainable development strategies and spatial land-use planning [16,17].
Traditional demographic statistics and analysis mainly rely on population surveys, including censuses and sampling studies. Until now, China has carried out six censuses. Although population surveys are scientific and authoritative [10,18], their data acquisition cycle is long and townships are the smallest survey unit, such that the spatial resolution of the data is insufficient [19]. Therefore most studies do not use the administrative unit as the research object [13,20,21]. With the rapid development of geographic information system and remote sensing technology, multi-source remote sensing data have been widely applied in spatial population research, especially land-use and night-time light (NTL) data [19,[22][23][24][25][26][27][28][29]. For example, Yang et al. [24] combined Defense Meteorological Satellite Program's Operational Linescan System (DMSP-OLS) NTL data, enhanced vegetation index data and digital elevation model (DEM) data to simulate the population density of Zhejiang. Hu et al. [25] determined the spatial distribution of population in Sichuan and Chongqing based on NTL data and land-use data. Other studies have shown that the spatial distribution of regional populations can be well-described by data processing, multi-source data fusion and model improvement [24,30]. Most research has remained focused on the spatial modeling of population at a single point in time [19,23,29,31,32], however few studies adopt multivariate data to model the population spatial distribution in a long time series.
Although the DMSP-OLS dataset provides continuous NTL data from 1992 to 2013 [33], its imagery contains problems due to OLS limitations such as discontinuity and oversaturation by bright lights [34,35] These data were replaced by Suomi National Polar-orbiting Partnership Visible Infrared Imaging Radiometer Suite (NPP-VIIRS) NTL data after 2013, bringing clear upgrades such as improved spatial resolution and reduced saturation [33] as well as on-board calibration [36]. Although these are clear upgrades, they also present challenges to obtaining consistent long-term NTL data [35,37,38], such that proper integration of the two datasets must be accomplished before the construction of a long-term population spatial distribution dataset and there have been several attempts to integrate DMSP and VIIRS NTL data [37][38][39][40][41]. Zhu et al. [39] established the relationship between the two at the provincial level and used it to model China's Gross Domestic Product. Zhao et al. [37] achieved this at the pixel level and established a long-term NTL dataset in Southeast Asia. Previous studies have contributed to enhancing the consistency of NTL between DMSP and VIIRS data, however there are limitations regarding a widespread application of current methods, such as the models proposed has regional limitations and may not be suitable for other regions [37,41]; the datasets used are not accessible to general public [40,41]; and the time series of data generated only has consistent NTL indices at the administrative level and is still limited at the pixel level [39].
The municipality of Chongqing integrates a metropolis and a large rural area that is mostly mountainous area, which is characterized by intense human activities and a fragile ecological environment. According to the NBSC, ongoing urbanization in Chongqing resulted in the rural population declining from 15.33 million to 10.7 million from 2005 to 2018 (a decrease of 30.2%), exceeding the national average of 24.34%. Meanwhile, Chinese policies targeting poverty alleviation and rural revitalization have benefitted most residents in poor mountainous areas through relocation, resulting in major changes in population distribution. Therefore, exploring the spatiotemporal changes in Chongqing's population via a timely understanding of population distribution data can help guide population migration from mountainous areas, promote the sustainable development of the regional economy and inform ecological restoration in mountainous areas.
This study explored integration methods for the two NTL datasets that are suitable for the study area at the pixel level and constructed a long-term NTL dataset that provides a basis for modeling long-term population spatial distribution data. Then it simulated the spatial distribution of Chongqing's population in 2000, 2005, 2010, 2015 and 2018 by integrating the two NTL datasets and analyzing spatiotemporal changes. Our results can serve as a scientific reference for rationally allocating urban and rural resources, optimizing urban and rural spatial patterns and promoting the high-quality development of the regional economy.

Study Area
Chongqing is located in the eastern Sichuan Basin, covering 8.24 × 10 4 km 2 from 28 • 10 -32 • 13 N and 105 • 11 -110 • 11 E (Figure 1). Its 26 districts and 12 counties cover a rugged landscape that is 75.33% mountainous. Its location at the intersection of the Silk Road and the Yangtze River Economic Belt allows it to form connections between east and west while driving economic development between north and south, leading to a vital role in China's development strategy underlain by the Belt and Road Initiatives and Yangtze River Economic Belt [42]. Chongqing is one of the important population areas in Western China, having a resident population of 31.02 million in 2018, an increase of 2.53 million compared with 2000 (NBSC); its average population density is about three times the national average. Its location in the upper reaches of the Yangtze River is part of an ecological protective screen within the Yangtze River Economic Belt. Its complex topography and fragile ecological environment enhance tensions between humans and the environment, such that ecological construction and regional development face many challenges. 2018 (a decrease of 30.2%), exceeding the national average of 24.34%. Meanwhile, Chinese policies targeting poverty alleviation and rural revitalization have benefitted most residents in poor mountainous areas through relocation, resulting in major changes in population distribution. Therefore, exploring the spatiotemporal changes in Chongqing's population via a timely understanding of population distribution data can help guide population migration from mountainous areas, promote the sustainable development of the regional economy and inform ecological restoration in mountainous areas. This study explored integration methods for the two NTL datasets that are suitable for the study area at the pixel level and constructed a long-term NTL dataset that provides a basis for modeling long-term population spatial distribution data. Then it simulated the spatial distribution of Chongqing's population in 2000, 2005, 2010, 2015 and 2018 by integrating the two NTL datasets and analyzing spatiotemporal changes. Our results can serve as a scientific reference for rationally allocating urban and rural resources, optimizing urban and rural spatial patterns and promoting the high-quality development of the regional economy.

Study Area
Chongqing is located in the eastern Sichuan Basin, covering 8.24 × 10 4 km 2 from 28°10′-32°13′N and 105°11′-110°11′E ( Figure 1). Its 26 districts and 12 counties cover a rugged landscape that is 75.33% mountainous. Its location at the intersection of the Silk Road and the Yangtze River Economic Belt allows it to form connections between east and west while driving economic development between north and south, leading to a vital role in China's development strategy underlain by the Belt and Road Initiatives and Yangtze River Economic Belt [42]. Chongqing is one of the important population areas in Western China, having a resident population of 31.02 million in 2018, an increase of 2.53 million compared with 2000 (NBSC); its average population density is about three times the national average. Its location in the upper reaches of the Yangtze River is part of an ecological protective screen within the Yangtze River Economic Belt. Its complex topography and fragile ecological environment enhance tensions between humans and the environment, such that ecological construction and regional development face many challenges.

Data Sources
DMSP-OLS Version 4 NTL data from 2000 to 2013 were obtained from the Paynes Institute for Public Policy, Colorado School of Mines (https://eogdata.mines.edu/dmsp/ downloadV4composites.html). These have a spatial resolution of 30 arc-seconds, with data values ranging from to 0-63 and have been denoised [43]. Monthly VIIRS Cloud Mask (vcm) data from 2013 to 2018 were also obtained from the Paynes Institute for Public Policy, Colorado School of Mines (https://eogdata.mines.edu/dmsp/download_radcal.html), with a spatial resolution of 15 arc-seconds that excludes observations affected by stray light. The data contained additional noise from sources such as auroras, fires, boats, other temporary lights and outliers, probably caused by stable lights from oil or gas fires.
Land-use data (1 km × 1 km) were obtained from the Resource and Environment Science Data Centre of the Chinese Academy of Sciences (http://www.resdc.cn/) with major categories including cultivated land, forest, grassland, water, residential land and unused land. Resident population data at the county level were obtained from the Chongqing Statistical Information Net (http://data.tjj.cq.gov.cn/), while those at the township level in 2015 were derived from the China County Statistical Yearbook 2016. DEM data were obtained from the Geospatial Data Cloud (http://www.gscloud.cn/).

Data Preprocessing
Firstly, all data were extracted by administrative boundaries and the DMSP-OLS NTL data was resampled to 1 km grids. Secondly, a stepwise calibration approach at the global scale was used to improve the temporal inconsistency of DMSP time series [44]. Thirdly, calculated the average value of VIIRS data from January to December to generate annual time series of VIIRS NTL imagery. Fourthly, mask extraction was then used to remove noise from the NPP-VIIRS NTL data. The DMSP-OLS NTL data and the annual NTL data provided by the NPP-VIIRS dataset were used as mask data. Masks were selected for each year according to the principle of time adjacency [34,45,46]. Finally, The maximum value of VIIRS NTL data in the main urban area of Chongqing was selected as the effective light intensity threshold and the eight-neighborhood algorithm was used to smooth the VIIRS NTL data [47]. These procedures allowed the NTL correction data to be obtained.

Methods
We established a relationship model between the two kinds of NTL data (based on the pixel scale), constructed a long time series of stable NTL datasets, then modeled the spatiotemporal dynamics of Chongqing's population from 2000 to 2018.

Integrating DMSP-OLS and NPP-VIIRS NTL Data
In order to match the spatial resolution and radiation characteristics of the two NTL data, we first performed two processes on the VIIRS data with reference to Zhao et al. [37]. One is using a kernel density (KD) method for resampling to make the spatial resolution the same as the DMSP data. The other is the logarithmic transformation. On this basis, we further discuss the NTL integration model and convert the value of VIIRS data.
(1) Spatial Resampling Using a KD Method Given that the blur of DMSP NTL image is a Gaussian point-spread function, the influence of neighborhood NTL brightness should be taken into account during the conversion of VIIRS spatial resolution. This paper adopted a quartic kernel function to realize as follows: where f(x) denotes the estimation of the KD function; n is the total number of samples; h is the window width and the value is five times of VIIRS pixel size here; K is the KD function; X is the pixel to be corrected; and X i is the neighbor pixels within the window.
(2) Logarithmic Transformation Logarithmic transformation of NPP-VIIRS data can better suppress the sharp radiance jump within urban core areas and strengthen the radiance variance within suburban and rural areas [37]. Therefore, we performed a logarithmic transformation for VIIRS data as follows: where N i denotes s the aggregation results of VIIRS NTLs using the KD method and Log_N i denotes the corresponding logarithmic transformation results. To avoid invalid values caused by logarithmic transformation, a constant of 1 was added. Considering that a slight seasonal difference may exist in annual VIIRS data, 2013 data were used to determine the relationship between the two data sets. We observed a positive correlation between DMSP and processed VIIRS value ( Figure 2).
where f(x) denotes the estimation of the KD function; n is the total number of samples; h is the window width and the value is five times of VIIRS pixel size here; K is the KD function; X is the pixel to be corrected; and Xi is the neighbor pixels within the window.
(2) Logarithmic Transformation Logarithmic transformation of NPP-VIIRS data can better suppress the sharp radiance jump within urban core areas and strengthen the radiance variance within suburban and rural areas [37]. Therefore, we performed a logarithmic transformation for VIIRS data as follows: where Ni denotes s the aggregation results of VIIRS NTLs using the KD method and Log_Ni denotes the corresponding logarithmic transformation results. To avoid invalid values caused by logarithmic transformation, a constant of 1 was added.
(3) Conversion of the VIIRS NTL Value Both DMSP and VIIRS products provide NTL data in 2012 and 2013 and the monthly VIIRS data in 2013 include all months, while the monthly data in 2012 are only available from April to December. Considering that a slight seasonal difference may exist in annual VIIRS data, 2013 data were used to determine the relationship between the two data sets. We observed a positive correlation between DMSP and processed VIIRS value ( Figure 2).
For further analysis, we developed a linear regression model, a quadratic polynomial regression model and a power function regression model relating the DMSP and processed VIIRS values in 2013, in order to find the best model for integrating NTL data.

Modeling the Spatiotemporal Dynamics of Population
NTL mainly comes from household lighting, roads, urban lightscapes, all of which are closely related to human activities. Moreover, NTL intensity directly reflects the intensity of such activities. Figure 3 shows the relationships between population density and the mean value of NTL at the county level. In Chongqing, NTL intensity grew rapidly with growth population growth. The quadratic polynomial model had the highest coefficient of determination of all models tested including the linear model and the power function model. For further analysis, we developed a linear regression model, a quadratic polynomial regression model and a power function regression model relating the DMSP and processed VIIRS values in 2013, in order to find the best model for integrating NTL data.

Modeling the Spatiotemporal Dynamics of Population
NTL mainly comes from household lighting, roads, urban lightscapes, all of which are closely related to human activities. Moreover, NTL intensity directly reflects the intensity of such activities. Figure 3 shows the relationships between population density and the mean value of NTL at the county level. In Chongqing, NTL intensity grew rapidly with growth population growth. The quadratic polynomial model had the highest coefficient of determination of all models tested including the linear model and the power function model. Different land-use patterns reflect population distribution and human production [48]. Our correlation analysis between population and land-use types at the county level showed that population was positively correlated with cultivated land, water and residential land at the 1% significance level, positively correlated with unused land at the 5% significance level and negatively correlated with forest and grassland at the 0.05 significance level (  Different land-use patterns reflect population distribution and human production [48]. Our correlation analysis between population and land-use types at the county level showed that population was positively correlated with cultivated land, water and residential land at the 1% significance level, positively correlated with unused land at the 5% significance level and negatively correlated with forest and grassland at the 0.05 significance level (Table 1). Table 1. Correlation analysis between population and land-use type (by area).
The population spatial distribution pattern was therefore closely related to NTL and land use, so we used the random-effect model to establish the relationship between population, NTL and land use. The resident population in each district and county was selected as the dependent variable and the total value of NTL (NT), the number of bright pixels (NL) and the number of dark pixels (ND) of each land use type in each district and county were used as independent variables. Considering that geographical factors are also important factors affecting population distribution, elevation variables were also added into the model as independent variables, which includes the number of pixels with altitudes (NPA) of 0-300 m, 300-500 m, 500-1000 m and >1000 m in each district and county. Prior to empirical simulation, stepwise regression was used to identify the key independent variables with a significance level within 20%. The key independent variables included the NT of cultivated land and forest; NL of residential; ND of cultivated land and grassland; and the NPA of 0-300 m, 300-500 m, 500-1000 m and > 1000 m. The collinearity between the variables was then tested using the variance inflation factor (VIF); the maximum VIF of a single variable was 3.42 and the overall VIF was 2.49 and they were well below the critical value of 10, indicating no serious collinearity problem between the variables. The empirical model settings were as follows: where Pit is the resident population of the ith county in the tth year; i = 1,2,…,38; t represents the known year; xit represents the observation value of variables in the ith county in the tth year; αi is the individual difference between regions; βi is a parameter to be  The population spatial distribution pattern was therefore closely related to NTL and land use, so we used the random-effect model to establish the relationship between population, NTL and land use. The resident population in each district and county was selected as the dependent variable and the total value of NTL (NT), the number of bright pixels (NL) and the number of dark pixels (ND) of each land use type in each district and county were used as independent variables. Considering that geographical factors are also important factors affecting population distribution, elevation variables were also added into the model as independent variables, which includes the number of pixels with altitudes (NPA) of 0-300 m, 300-500 m, 500-1000 m and >1000 m in each district and county. Prior to empirical simulation, stepwise regression was used to identify the key independent variables with a significance level within 20%. The key independent variables included the NT of cultivated land and forest; NL of residential; ND of cultivated land and grassland; and the NPA of 0-300 m, 300-500 m, 500-1000 m and >1000 m. The collinearity between the variables was then tested using the variance inflation factor (VIF); the maximum VIF of a single variable was 3.42 and the overall VIF was 2.49 and they were well below the critical value of 10, indicating no serious collinearity problem between the variables. The empirical model settings were as follows: where P it is the resident population of the ith county in the tth year; i = 1,2, . . . , 38; t represents the known year; x it represents the observation value of variables in the ith county in the tth year; α i is the individual difference between regions; β i is a parameter to be estimated; and µ it is a random error term. Considering the real situation of population distribution, water and unused land were not involved in the model calculation [19]. Next, based on the estimated results of the random-effect model, the resident population of each grid was calculated as follows: where P ijk is the resident population in the kth pixel of the jth land use type in the ith county; P 0 is a constant; N i is the number of pixels in the ith county; a j , b j , c j , and d j are coefficients; m is the number of land-use types; NT ijk , NL ijk and ND ijk are the total value of NTL, the number of bright pixels and the number of dark pixels in the kth pixel of the jth land use type in the ith county, respectively; NPA ikn is the number of pixels of the nth elevation interval in the kth pixel and ith district and county. Negative coefficients for some variables in the simulation equation established by the random-effect model resulted in the estimated population of some pixels being negative, a situational impossibility. Therefore, pixels with a negative estimation value were assigned a value of 0 before obtaining the preliminary estimated population data. Finally, the statistical data for county population were used to adjust the simulation results as follows: where P ijk is the final resident population in the kth pixel of the jth land use type in the ith county; P i is the statistical data of the resident population in the ith county; and P i is the total population by preliminary estimate in the ith county.

Evaluation of Model Accuracy
Based on the population census data at the township level, the correlation coefficient (R), mean absolute error (MEA), mean relative error (MRE) and root mean square error (RMSE) were selected to evaluate accuracy as follows: where P i is the statistical resident population in the ith township provided by census data, PE i is the estimated resident population in the ith township, P is the average of the statistical population and PE and is the average of the estimated population.

Integration Model
The power function model had the highest coefficient of determination (R 2 = 0.907) of the three models tested (Figure 4). Therefore, the relationship established by the power function model was used to simulate DMSP data from 2014 to 2018. The method of integrating NTL data is as follows: TNL n = TNL a where TNL a n is the NTL radiance value for the DMSP-OLS data in the nth year; Log − N n is the processed VIIRS radiance value in the nth year; and TNL n is the value for the NTL integration data in the nth year.
where is the NTL radiance value for the DMSP-OLS data in the nth year; − is the processed VIIRS radiance value in the nth year; and is the value for the NTL integration data in the nth year.

Accuracy Assessment
We assessed the accuracy of the integrated NTL data by comparing the DMSP-OLS and adjusted NPP-VIIRS data in 2013 ( Table 2). The MRE value of the mean NTL generated from the integrated data was 8.19%, while the relative error (RE) values of the mean NTL varied by county, with 39.47% of counties underestimated, 60.53% overestimated and 71.05% having RE values within 10%. The maximum and minimum RE values were 42.46% (Chengkou) and 0.2% (Yubei).
A previous study exploring the relationship between the two kinds of NTL data at the provincial level produced an MRE of 12.97% [39]. In comparison, our methods clearly improved the matching accuracy, making this approach feasible for integrating NTL data.

Accuracy Assessment
We assessed the accuracy of the integrated NTL data by comparing the DMSP-OLS and adjusted NPP-VIIRS data in 2013 ( Table 2). The MRE value of the mean NTL generated from the integrated data was 8.19%, while the relative error (RE) values of the mean NTL varied by county, with 39.47% of counties underestimated, 60.53% overestimated and 71.05% having RE values within 10%. The maximum and minimum RE values were 42.46% (Chengkou) and 0.2% (Yubei). A previous study exploring the relationship between the two kinds of NTL data at the provincial level produced an MRE of 12.97% [39]. In comparison, our methods clearly improved the matching accuracy, making this approach feasible for integrating NTL data.

Random-Effect Model
The random-effect model results produced an overall F value of 224.42, an R 2 value between groups of 0.71 and an overall p-value of 0.000, indicating that the model was well-established and that the modeling equation was reasonable (Table 3). Note: (1) *, ** and *** are significantly different from zero at the 10%, 5% and 1% levels, respectively.

Accuracy Assessment
We evaluated the population modeling results using 2015 census data for 150 randomly selected villages and towns. As terrain factors could affect the accuracy of population simulations, we divided the study area into the three zones by elevation (high-altitude, ≥1000 m; medium-altitude, 500-1000 m; and low-altitude, <500 m) among which the randomly selected villages and towns were evenly distributed ( Figure 5).

Spatial Distribution of Population in Chongqing
According to the fifth census in 2000, if the population density of a municipal district was more than 1500 persons/km 2 , the entire population was classified as urban. On this basis, we regarded population densities of > 1500 persons/km 2 as high-population-density regions. In addition, according to Tan et al.'s [19] hierarchical classification method for population density, areas with a population density of 200-1500 people/km 2 were classified as intermediate-density regions and areas with a population density < 200 people/km 2 were classified as low-density regions.
From 2000 to 2018, Chongqing's population density has generally increased in the west and decreased in the east (Figure 6). High-density regions were mainly distributed in western Chongqing and those centered on Yuzhong continued to expand. In contrast, the population density of most regions in the northeast and southeast decreased to varying degrees, trending toward low population density.

Spatial Distribution of Population in Chongqing
According to the fifth census in 2000, if the population density of a municipal district was more than 1500 persons/km 2 , the entire population was classified as urban. On this basis, we regarded population densities of >1500 persons/km 2 as high-population-density regions. In addition, according to Tan et al.'s [19] hierarchical classification method for population density, areas with a population density of 200-1500 people/km 2 were classified as intermediate-density regions and areas with a population density <200 people/km 2 were classified as low-density regions.
From 2000 to 2018, Chongqing's population density has generally increased in the west and decreased in the east (Figure 6). High-density regions were mainly distributed in western Chongqing and those centered on Yuzhong continued to expand. In contrast, the population density of most regions in the northeast and southeast decreased to varying degrees, trending toward low population density.
Low-density regions in Chongqing grew from 35.64 × 10 3 km 2 in 2000 to 41.07 × 10 3 km 2 in 2018 (an increase of 15.22%) ( Table 5). 95% of the newly added regions were created by the loss of population from Intermediate-density regions, mainly in the northeast and southeast, including Fengjie, Yunyang, Wushan, Wuxi, Xiushan, Fengdu and Shizhu (Figure 7).   Table 5). 95% of the newly added regions were created by the loss of population from Intermediate-density regions, mainly in the northeast and southeast, including Fengjie, Yunyang, Wushan, Wuxi, Xiushan, Fengdu and Shizhu (Figure 7). Table 5. Changes in population density from 2000 to 2018.  The total intermediate-density area decreased from 45.02 × 10 3 km 2 in 2000 to 38.38 × 10 3 km 2 in 2018 (a decrease of 14.75%). The reduced regions were mainly distributed in the primary urban zone of Chongqing and in the northeast. Intermediate-density regions in urban zone tended to agglomerate and gradually develop into high-density regions, while intermediate-density regions in the northeast gradually lost their population and developed into low-density regions. In addition, within each district and county, population development trended toward agglomeration, manifested as a gradual increase in urban population density and the gradual evolution of intermediate-density regions into highdensity regions; however, in non-urban areas, population loss was more common in intermediate-density regions, where population density decreased.  The total intermediate-density area decreased from 45.02 × 10 3 km 2 in 2000 to 38.38 × 10 3 km 2 in 2018 (a decrease of 14.75%). The reduced regions were mainly distributed in the primary urban zone of Chongqing and in the northeast. Intermediate-density regions in urban zone tended to agglomerate and gradually develop into high-density regions, while intermediate-density regions in the northeast gradually lost their population and developed into low-density regions. In addition, within each district and county, population development trended toward agglomeration, manifested as a gradual increase in urban population density and the gradual evolution of intermediate-density regions into high-density regions; however, in non-urban areas, population loss was more common in intermediate-density regions, where population density decreased.

Regional Division
High-density regions gradually expanded from 1.07 × 10 3 km 2 in 2000 to 2.28 × 10 3 km 2 in 2018 (an increase of 113.08%). In 2000, these were mainly distributed within a radius of 24 km from Yuzhong ( Figure 8) but ongoing urbanization expanded this range to a radius of 33 km by 2018. In addition, urban areas within each district and county also became distributed within high-density regions, which expanded to different degrees. The low-altitude zone had the highest average population density and population growth while trending toward agglomeration (Table 6). The average population density here increased from 550.58 to 647.08 people/km 2 from 2000 to 2018, a total population increase of 3.16 million. The growth rate was fastest from 2010 to 2015, increasing by 134.69 × 10 4 people in only five years. In contrast, the medium-and high-altitude zones showed declining population density and total population. The medium-altitude zone showed a drop in average population density from 232.66 to 223.70 people/km 2 from 2000 to 2018, a total population decrease of 0.4 million. The average population density in the high-altitude zone dropped from 64.45 to 58.78 people/km 2 from 2000 to 2018, a total population decrease of 0.18 million. Table 6. Average population density and total population of each altitude zone.

Discussion
The DMSP-OLS dataset represents the most widely used NTL data over the previous two decades, while the new NPP-VIIRS NTL data have been available since 2012. Despite the great significance of studying long-term population evolution in the context of urbanrural migrations, few studies have integrated the two datasets to simulate and monitor population spatial changes over the full time period. In this study, we proposed a method for integrating the DMSP-OLS and NPP-VIIRS data at the pixel scale in order to extend the temporal coverage of NTL data. Meanwhile, we have evaluated the accuracy of the integrated NTL data and the MRE was 8.19%. Our integration accuracy was improved by 4.78% compared with the long-time-series NTL dataset established at the provincial level [39], which indicated that our method for NTL integration was feasible and the resulting data had good quality and generally reliable temporal consistency.
Previous studies have simulated population spatial distribution in different regions using NTL and land-use data. Hu et al. [25] did this for Sichuan and Chongqing in 2014, with MREs for population data based on DMSP-OLS and NPP-VIIRS NTL data of 46.3% and 44.62%, respectively. Chowdhury et al. [23] developed a model for estimating the population in the Indian portion of the Indo-Gangetic Plains at both city and state levels by employing OLS NTL data. The model was validated for the population of year 1995, with an MRE of 9.4%. Liu et al. [26] simulated the spatial pattern of urban and rural residents in the Huang-Huai-Hai area with an MRE of 15.6%. Tan et al. [19] simulated the population density of China in 2000, achieving a correlation coefficient between the statistical and simulated values of 0.95. The accuracy of population simulations in mountainous areas such as Chongqing and Sichuan is lower than in plains areas such as Huang-Huai-Hai, demonstrating that population simulation in mountainous areas is more challenging and uncertain. As we were limited by the difficulty of obtaining accurate population data in towns and villages, we only tested the accuracy of population simulation in 2015; the R value (0.85) and MRE (26.98%) confirmed that the adjusted VIIRS data were capable of effectively simulating spatial population patterns. We optimized the simulation method for mountainous areas based on previous research [25], increasing the results' accuracy by nearly 20%. We also introduced a feasible method for constructing long-term population spatial data, which is helpful for scientifically monitoring spatiotemporal trends in mountainous populations. In addition, the U.S. Department of Defense has developed the Landscan database using an innovative approach with Geographic Information Systems and Remote Sensing, which is the finest resolution global population distribution data available [49]. In order to further verify our results, we also evaluated Landscan data using 2015 census data for 150 randomly selected villages and towns and the results showed that the R value and MRE were 0.78 and 35.7% respectively, which also proved the feasibility of our method.
It is worth mentioning that there are still some limitations in this study. First, although we were able to improve the accuracy of mountainous population spatial simulation through data processing, this method was unable to completely eliminate inherent defects in the DMSP-OLS data, such as light saturation in urban centers with high light intensity [50] and insufficient detection capabilities in low-radiation areas such as rural areas [33]. These flaws reduce the accuracy of population simulation to a certain extent. Second, the change of lighting technology (from sodium vapor to light-emitting diode) reduced NTL values in the city center [51], which may have led to an underestimation of population simulation results. Third, the study was difficult to obtain the annual population distribution data and we only simulated the population distribution in the five periods of 2000, 2005, 2015 and 2018 due to limitation of data collection. Fourth, compared with DMSP-OLS data, NPP-VIIRS data have a higher spatiotemporal resolution. The advantages of the latter were not fully integrated into the long-term NTL dataset and further research is needed to improve the spatial resolution of NTL integration.

Conclusions
We integrated DMSP-OLS and NPP-VIIRS NTL data to construct a long-term NTL dataset, using the random-effect model with land-use data and corrected NTL data to model the spatiotemporal dynamics of the Chongqing's population from 2000-2018. At the pixel level, there was a power function relationship between the two datasets (R 2 = 0.907). Compared with an NTL integration model previously established at the provincial level, our model was 4.78% more accurate. In addition, accuracy tests using 2015 data resulted in an MRE of 26.98%, an improvement of nearly 20% when compared with previous studies of mountainous populations. Therefore, our approach is feasible and provides a technical method for monitoring spatiotemporal population changes in mountainous areas.
From 2000-2018, the spatial distribution of Chongqing's population has increased in the west and decreased in the east, while also increasing in low-altitude areas and decreasing in the medium-high altitude areas. Moreover, population agglomeration was common. At the provincial level, high-density regions showed a significant increase, while decreasing in intermediate-density regions. The population density significantly increased in the central urban area and immediate surroundings in every district and county, while significantly decreased in non-urban areas, especially in the northeast.