Optimization of Modelling Population Density Estimation Based on Impervious Surface

Population data are key indicators of policymaking, public health, and land use in urban and ecological systems; however, traditional censuses are time-consuming, expensive, and laborious. This study proposes a method of modelling population density estimations based on remote sensing data in Hefei. Four models with impervious surface (IS), night light (NTL), and point of interest (POI) data as independent variables are constructed at the township scale, and the optimal model was applied to pixels to obtain a finer population density distribution. The results show that: (1) impervious surface (IS) data can be effectively extracted by the linear spectral mixture analysis (LSMA) method; (2) there is a high potential of the multi-variable model to estimate the population density, with an adjusted R2 of 0.832, and mean absolute error (MAE) of 0.420 from 10-fold cross validation recorded; (3) downscaling the predicted population density from the township scale to pixels using the multi-variable stepwise regression model achieves a more refined population density distribution. This study provides a promising method for the rapid and effective prediction of population data in interval years, and data support for urban planning and population management.


Introduction
Population data are considered important indices for the development of a country or region. Urban or land planning, policymaking, public emergencies, and other public aspects require detailed population data [1,2]. Currently, most countries obtain population data through the census. However, accessing the population data is difficult because of the dispersion and dynamics of the population distribution [3,4]. In China, the census is the most common method for collecting population data. Although the census is relatively credible, it is time-consuming, costly, and obstructed by difficulties in reality [5]. Nowadays, population mobility makes it difficult to obtain actual residence figures. People are more protective of their privacy. It is very common to refuse to fill in personal information. The situations above increase the difficulty of the census. Moreover, the national census is conducted once every ten years, which may lead to the absence of population data in interval years. Consequently, scientific decision making lacks support.
Unlike the national census, remote sensing data are not subject to time restrictions [6][7][8][9]. Owing to the rapid development of remote sensing and computer technologies, some presented study is to improve the population density estimation based on impervious surfaces. To better delineate the population distribution area, our study integrates POI data and NTL data with IS data to improve the model. Therefore, in this study, the geo-spatiotemporal big data POI and the latest highresolution NTL data from Luojia-1 were fused with the IS information extracted from remote sensing images to try to establish the population estimation model. With IS data as the main independent variable, the population distribution of Hefei was explored using the stepwise regression method. The optimal population estimation model was then applied to the pixel scale. Finally, a finer population distribution map was obtained.

Study Area
Hefei (30 • 57 -32 • 32 N, 116 • 41 -117 • 58 E; Figure 1) is situated in the middle of Anhui Province in China, covering a total area of 11,408.48 km 2 with an average altitude of 30 m. The terrain comprises plains and low hills. It is characterised by a subtropical monsoon climate with an annual mean daily temperature of 16 • C and a total annual precipitation of 995.2 mm. As the capital, Hefei is the cultural, commercial, financial, and political centre of Anhui Province. Located in the radiation belt of the Yangtze River Delta, it acts as a gateway for the development of the central and western regions. Meanwhile, as a node city of the 'One Belt, One Road' strategy, Hefei has great economic potential [41].

Data Collection
The main data used in this study include the satellite data, POI data, NTL data, township vector boundary data in Hefei, and census data, for 2018 (Table 1). Landsat 8 OLS multi-temporal images were selected to extract IS data on April 10, 2018 because of the sparse vegetation in spring but lush plants in summer and a large area of bare soil in winter affecting the experiment. The cloudiness of this image was 0.03, which had negligible impacts on data processing. Radiometric calibration, atmospheric correction, and geo-referencing were conducted from the remote sensing image downloaded for free from the United States Geological Survey (USGS, http://earthexplorer.usgs.gov/, accessed on 15 April 2020). Subsequently, the study area was masked from the image using administra- stepwise regression model for population estimation was applied to provide a reference for mapping the population density distribution.

Data Collection
The main data used in this study include the satellite data, POI data, NTL data, township vector boundary data in Hefei, and census data, for 2018 (Table 1). Landsat 8 OLS multi-temporal images were selected to extract IS data on April 10, 2018 because of the sparse vegetation in spring but lush plants in summer and a large area of bare soil in winter affecting the experiment. The cloudiness of this image was 0.03, which had negligible impacts on data processing. Radiometric calibration, atmospheric correction, and geo-referencing were conducted from the remote sensing image downloaded for free from the United States Geological Survey (USGS, http://earthexplorer.usgs.gov/, accessed on 15 April 2020). Subsequently, the study area was masked from the image using administrative data from the Hefei administrative boundary. Census data were obtained from the China County Statistical Yearbook published in 2019, and the Hefei Municipal Bureau Statistics. The 141 townships were used for modelling to explore the population density estimation model at a small scale. Moreover, the POI data, including 123,348 POIs in 2018, were crawled from the Baidu Map Services (http://map.baidu.com, accessed on 3 January 2021), which is the most widely used and largest web map service provider in China [7].
The NTL data of Luojia-1 were downloaded from the High-Resolution Earth Observation System of Hubei Data and Applications Network (http://www.hbeos.org.cn/, accessed on 27 January 2021). Finally, all spatial data were georeferenced in the same projection owing to the different source data. The flowchart of modelling the population density estimation is shown in Figure 2. (http://map.baidu.com, accessed on 3 January 2021), which is the most widely used and largest web map service provider in China [7]. The NTL data of Luojia-1 were downloaded from the High-Resolution Earth Observation System of Hubei Data and Applications Network (http://www.hbeos.org.cn/, accessed on 27 January 2021). Finally, all spatial data were georeferenced in the same projection owing to the different source data. The flowchart of modelling the population density estimation is shown in Figure 2.

Mapping IS Distribution and Validating the Results
LSMA can effectively discriminate different objects, such as soil, water, vegetation, and ISs. Therefore, LSMA was utilised in the mapping process to obtain the proportion of an IS within a pixel.

Mapping is Distribution and Validating the Results
LSMA can effectively discriminate different objects, such as soil, water, vegetation, and ISs. Therefore, LSMA was utilised in the mapping process to obtain the proportion of an IS within a pixel.
After pre-processing, the image was separated from other land covers using the modified normalised difference water index (MNDWI) [42], aiming to remove the non-IS fractions first. Subsequently, the minimum noise fraction (MNF) was adopted to reduce the redundancy within the image and improve the quality of the endmember selection [29]. Generally, three or four endmembers are selected [43]. As a consequence of confusing pixels in ISs, we selected four endmembers through repeated trials, including high-and low-albedo objects, vegetation, and bare soil, based on the vegetation-impervious surfacesoil (V-I-S) model [44]. Finally, the fully constrained LSMA method was used to develop the fractional IS map, during which the high-and low-albedo objects were added together [45]. Thus, a map of the IS distribution in 2018 was obtained.
When measuring the V-I-S model fitness, the average root mean square (RMS) of overall bands was used to evaluate the accuracy of the LSMA (Equation (1)) [44]: where RMS is the root mean square, m is the number of pixels in this image and ε(λ i ) is the residual error of the pixels after unmixing.
To further assess the unmixing accuracy, detailed ground reference data (the actual IS data) are required for comparison with the experimental data [46]. In this study, a total of 120 validated samples were randomly selected from satellite images in Hefei, each containing 3 × 3 pixels (90 × 90 m). Subsequently, all the validation samples were overlapped on the high-resolution image in 2018 from Google Earth Pro, which has been proven suitable for obtaining ground reference data [46]. Finally, the linear fit between the true and experimental data was evaluated using the coefficient of determination (R 2 ) and the root mean square error (RMSE). As the study is intended to estimate the population density using the IS data as a variable, the proportion of an IS in a pixel was calculated, and the mean data were used as a proxy for each township.

Data Preparation for Modelling
At present, the commonly used scales are 250, 500, and 1000 m in the study of population spatialisation. The selection of scale is mainly based on the area of the target and model stability. The 250 m scale and below is suitable for research in villages and communities, 500 m for that in counties and cities, and 1000 m for that in large-scale regional studies involving cities and provinces. According to existing research, the number of samples exceeding 50,000 affects the stability of the regression model and the accuracy of the model [13,40]. Therefore, the 500m scale was selected for modelling.
In this study, each 500 m grid was defined as a pixel. The pixel was built using ArcGIS 10.3 software. Using the 'Create Fishnet' function, the grid files of 500 m pixels were generated within the Hefei administrative boundary. The following experiments were conducted using ArcGIS 10.3 software.
The IS data were resampled to a resolution of 500 m. We obtained population-related POI data from Baidu Map in 2018 using Python to crawl the eight categories, including restaurants, shopping, hospitals and clinic facilities, education facilities, entertainment and retail, public service facilities, companies, and residential areas [7]. To avoid multicollinearity caused by multiple variables, the POI data of eight categories were merged into an image.
The commonly used point-based methods include the analysis based on Euclidean distance and point density. The Euclidean distance considers that the plane space is homogeneous, ignoring the relationship of facilities and service functions between cities, which is suitable for studying scattered points. Point density analysis is based on the aggregation of point distribution, and these points often interact with each other in space. We tried to aggregate the POIs into the grids to obtain the POI density. However, there was no distribution of POIs in the grids of sparse population areas, leading to the null data in samples.
The kernel density analysis is a point density method based on the first law of geography [7] and is suitable for the estimation of continuous geographic phenomena. POIs are distributed in the range of human activities and interact with urban facilities and transportation, and are consistent with aggregation distribution and continuous geographic phenomena.
Kernel density analysis was then conducted to analyse the point density spatially and discern the hotspots [47][48][49]. The crucial step was selecting the appropriate bandwidths [48]. We tested bandwidths varying from 1000 to 6000 m at an interval of 100 m, during which the correlation between the POI density generated by different bandwidths and the actual population density was analysed. Finally, it was found that when the bandwidth was 3500 m, the correlation between POI data and population density was the strongest. Therefore, the bandwidth of 3500 m performed well in identifying the hotspots, and we obtained the raster map of POI data at a 500 m pixel scale. The NTL image was clipped and pre-processed after downloading and subsequently resampled to pixels.
The vector file of the township boundary was superimposed with the raster image to obtain the township information on each pixel. To achieve consistency in the dimensions of different variables, the maximum and minimum standardisation methods were adopted to standardise the values of the three independent variables from 0 to 1 [21]. Logarithmic transformation was used to recalculate the population density data to eliminate the negative effects of the excessive population density [7]. Considering each township as a sample, the average value of each variable in each township was obtained as the independent variable value of a sample.

Model Concept and Validation
ISs are known to be closely related to population data [4]. However, the result of using IS data as the sole independent variable for population density estimations always leads to a coarse prediction, as a single variable cannot satisfactorily describe a populated zone [22,23]. Therefore, POI and NTL data were then added to the model. Assuming that the IS data has no significant multicollinearity with the other two data sources, the population density estimation model can be constructed (Equation (2)) [45]: where Y is the dependent variable (population density), X i is the independent variable (IS, POI, and NTL data), β 0 is the intercept, β i is the regression coefficient of the ith independent variable, ε is the error of the models, and n represents the number of independent variables. To explore the relationship between IS data, NTL data, and POI data, Pearson correlation and partial correlation analyses were conducted. The former aims to measure the intensity of the monotonic relationship between variables [50], and the latter can analyse the correlation between two variables by keeping the other variables constant and eliminating their effects [51]. The correlation analysis provided the basis for modelling. Furthermore, to quantitatively characterise the multicollinearity, it was checked according to a rule of thumb stating that a variance inflation factor (VIF) value above 10 rules out the variable because of the high multicollinearity [52].
The models were constructed by randomly selecting samples from the townships. Hence, a stratified sampling procedure was conducted. To better assess the accuracy of the model, the 10-fold cross validation was applied to evaluate the performance of all models [7]. Approximately 90% of the townships (divided into nine groups) were randomly selected from all 141 samples, and the rest (one group) were used to repeat 5 trials for the 10-fold cross validation. The census data of nine groups were used for training samples and one group for accessing the accuracy of the results. As the validated model can be applied for the population estimation, it is vital to test the model fit. In this study, the model fit between the predicted and true population density at the township scale was examined using the residual-based adjusted R 2 , relative mean square error (RMSE), and the mean average error (MAE) [52].
As the independent variable data used in the modelling at the township scale were the average independent variable values of all pixels in the township, the optimal model can directly be applied to the pixel [4]. The IS data of each pixel were extracted by LSMA and then resampled, NTL data were obtained by directly processing the Luojia-1 NTL data image, and the POI data were obtained by calculating the kernel density. The population density data and independent variables of each township were obtained by connecting the township name with the 'Identity' function of ArcGIS 10.3 software.
The estimation models are based on stepwise linear regression. Despite its limitations, stepwise regression is still very popular in recent studies [53][54][55], because it can effectively identify the best variables in many related or unrelated variables to build a good prediction model [56,57]. Finally, the IS, NTL, and POI data were merged into the pixels using their IDs for predicting the population density in a pixel with the coefficients of the optimal model. To validate the population density estimation model at the 500 m pixel scale, three administrative units were selected as proxies for high-, medium-, and low-density areas, respectively.

Analysis of is Distribution and Assessment
As the main independent variable of the model, the accuracy of the ISs directly affects the initial regression result. After calculating the MNDWI index, the water was masked. Four endmember types were distinguished based on the terrain properties, and finally, the high-and low-albedo endmembers were superimposed to obtain the distribution map of the IS proportion ( Figure 3). The central part of Figure 3 shows the four main districts of Hefei, wherein the deeper colour indicates a higher IS value. Almost every township has a deep-coloured area, representing a population gathering and large IS distribution. The RMS was tested at 0.094 within a reasonable scope [58].
The distribution of the validation points is shown in Figure 4. After calculating the proportion of each endmember type in the external square of each point, the proportion of the IS in the square was compared with the IS abundance value extracted by LSMA. Subsequently, the experimental and true data were compared by calculating the linear fitting between them. The linear regression was up to the standard, with R 2 = 0.788 and the RMSE = 0.129 ( Figure 5). Notably, a strong correlation was observed between experimental and true data.

Correlation Analysis between Variables
Pearson correlation analysis was conducted between independent variables, and the corresponding coefficients were calculated ( Table 2). The IS data were positively correlated with the NTL data and POI data, with correlation coefficients of 0.729 and 0.689, respectively, and with all being less than 0.9, which is the threshold for correlation analysis in population estimation [59]. The partial correlation coefficient of each variable (Table 3) is smaller than the correlation coefficient.
Four endmember types were distinguished based on the terrain properties, and finally the high-and low-albedo endmembers were superimposed to obtain the distribution ma of the IS proportion (Figure 3). The central part of Figure 3 shows the four main district of Hefei, wherein the deeper colour indicates a higher IS value. Almost every townshi has a deep-coloured area, representing a population gathering and large IS distribution The RMS was tested at 0.094 within a reasonable scope [58]. The distribution of the validation points is shown in Figure 4. After calculating th proportion of each endmember type in the external square of each point, the proportio of the IS in the square was compared with the IS abundance value extracted by LSMA Subsequently, the experimental and true data were compared by calculating the linea fitting between them. The linear regression was up to the standard, with R 2 = 0.788 an the RMSE = 0.129 ( Figure 5). Notably, a strong correlation was observed between exper mental and true data.

Correlation Analysis between Variables
Pearson correlation analysis was conducted between independent variables, and the corresponding coefficients were calculated ( Table 2). The IS data were positively correlated with the NTL data and POI data, with correlation coefficients of 0.729 and 0.689 respectively, and with all being less than 0.9, which is the threshold for correlation analysis in population estimation [59]. The partial correlation coefficient of each variable (Table 3) is smaller than the correlation coefficient. Table 2. Pearson correlation coefficients between all variables in the township scale (***-at the 1% level, **-at the 5% level).

IS Data
POI Data NTL Data IS data 1 0.729 *** 0.689 *** POI data 1 0.673 ** NTL data 1 Figure 5. Result of the linear fitting of reference and experimental data of the IS. Table 2. Pearson correlation coefficients between all variables in the township scale (***-at the 1% level, **-at the 5% level).

IS Data POI Data NTL Data
IS data 0.495 *** 0.391 *** POI data 0.495 *** 0.345 *** NTL data 0.391 *** 0.345 *** To further explore the correlation between the three independent variables, the variance inflation factor (VIF) values were calculated when the population density was the dependent variable [57] (Equation (3)): where R i is the correlation coefficient of the ith independent variables.
The VIF values were all below 3 (Table 4). Based on the above analysis, the correlation between independent variables was weak and the multicollinearity was low, which provided a basis for constructing the multi-variable regression model [57]. Note that the bivariate model (a) included IS data and NTL data, and the bivariate model (b) included IS data and POI data.

Comparison and Validation of Models
In this study, the error between the predicted and actual population densities was calculated by mean absolute error (MAE) (Equation (4)): where MAE is the relative error of the population density estimation, Pop a is the actual data of population density on a log scale, Pop b is the prediction data of the population density on a log scale, and n is the number of townships. The 10-fold cross validation results show that the single variable regression model for population density estimation was built with IS data as a single independent variable to test its prominence, which performs well, with R 2 = 0.689 and the RMSE = 0.922 (Table 5). Alternatively, the addition of POI data changed the model fit from 0.689 to 0.834 and the RMSE from 0.922 to 0.686, indicating the influence of POI data on the population density. The model fit performed best after adding all three independent variables in the regression, and the RMSE continued to decrease slightly.  Table 5 shows the results of the mean value of the adjusted R 2 , RMSE, and MAE. It indicates that the value of MAE decreased with the introduction of NTL data and POI data, and decreased most obviously after the introduction of POI data. In the bivariate models, model (b) with IS data and POI data as independent variables performed better. The MAE of the multi-variable model after the introduction of NTL data and POI data was the lowest, and the model fit was the best, which is regarded as the optimal model. In the optimal model, the group of minimum MAE was selected to explore the relationship between the predicted population density and the actual population density ( Figure 6). The model achieved a cross-validated MAE of 0.420 and an adjusted R 2 of 0.832. Before population spatialization, it is necessary to modify the population density between the township scale and the pixel scale to reduce the error caused by downscaling  Before population spatialization, it is necessary to modify the population density between the township scale and the pixel scale to reduce the error caused by downscaling [60] (Equation (5)): where Pop i,j is the modified population density of jth pixel of the ith township after exponential transformation, Pop b_ i,j , is the prediction data of population density of the jth pixel of the ith township after exponential transformation, Pop a_i is the actual data of population density of ith township, N is the number of townships, Area j is the area of the jth pixel. Figure 7a shows the population distribution maps based on census data, and Figure 8 shows those based on predicted data. It can be seen from Figure 7a that the population density within the administrative unit is homogeneous, and distinguishing between densely populated and sparsely populated areas is not possible; however, the township scale is the smallest in the census. Figure 7b shows the population distribution using a multi-variable regression model at a 500 m pixel scale. In the main urban areas with large population density gaps, the population distribution is no longer divided by administrative boundaries, and the value of each pixel represents the population density in the area. Therefore, this map more precisely reflects the population distribution within the administrative unit.
The current study areas for population estimation are concentrated in countries at the province or county scale [23,25], always leading to a coarse population distribution. Our study explored the population density estimation model in Hefei (a developing city), emphasising the importance of population data for some policies, such as the introduction of talent and household registration for migrant workers, which is closely concerned with the population. The model for population estimation can predict the population in the interval years of census data and then be applied to the pixel scale to redistribute the population. Owing to the lack of census data at the pixel scale, it was difficult to assess the accuracy of the population spatialisation quantitatively. Therefore, we compared the predicted data with high-resolution images and maps of administrative units at the township scale with three representative areas: high-, medium-, and low-density areas (Figures 8-10) [60]. The current study areas for population estimation are concentrated in countries at the province or county scale [23,25], always leading to a coarse population distribution. Our study explored the population density estimation model in Hefei (a developing city), emphasising the importance of population data for some policies, such as the introduction Located at the junction of the main city of Hefei and Feixi County, Shangpai covers an area of 121 km 2 , with a population density of 2228.29 people/km 2 . From the high-resolution image (Figure 9a), the central area of Shangpai Township is located in the northern part, and several dwelling units are sprinkled in the south-eastern area. Other areas are mainly forests or croplands. The population density of the Shangpai town, displayed in Figure 9b, does not consider the internal population distribution differences. The yellow pixels shown in Figure 9c correspond to the building area in Figure 9a, which represents the area where the population is concentrated. Xiage Township belongs to Chaohu City (a county-level city in Hefei), covering an area of 184 km 2 , with a population density of 306.27 people/km 2 . It is located on the northern bank of Chaohu Lake, with beautiful scenery. The township (Figure 10a) mainly consists of mountainous terrain with vast tracts of forests, and the residential areas are distributed in the southern part, as shown in Figure 10c.  Located at the junction of the main city of Hefei and Feixi County, Shangpai covers an area of 121 km 2 , with a population density of 2228.29 people/km 2 . From the high-resolution image (Figure 9a), the central area of Shangpai Township is located in the northern part, and several dwelling units are sprinkled in the south-eastern area. Other areas are mainly forests or croplands. The population density of the Shangpai town, displayed in Figure 9b, does not consider the internal population distribution differences. The yellow pixels shown in Figure 9c correspond to the building area in Figure 9a, which represents the area where the population is concentrated. Xiage Township belongs to Chaohu City (a county-level city in Hefei), covering an area of 184 km 2 , with a population density of 306.27 people/km 2 . It is located on the northern bank of Chaohu Lake, with beautiful scenery. The township (Figure 10a) mainly consists of mountainous terrain with vast tracts of forests, and the residential areas are distributed in the southern part, as shown in Figure 10c. Figures 8-10 depict the comparison between Google Earth high-resolution images, census data, and population density at a 500 m pixel scale. Notably, all population distributions at the township scale were homogenous. From the finer-scale maps of the three figures, the pixels at the administrative boundary of the townships are segmented. The value of independent variables in the segmented pixel is determined by the township with the pixels' largest area. However, the two segmented parts are involved in the calculation, resulting in the repeated calculation of the boundary part. This is one error source in the pixel-scale population distribution map.

Discussion
Testing on extracting the proportion of ISs by LSMA obtained a satisfying result, corresponding with the previous findings [29,61]. The map of the ISs showed that their distribution was related to the population (Figures 3 and 7). However, estimating the population with a sole variable may led to a coarse result. Thus, POI and NTL data were considered as added independent variables. Note that the number of independent variables in the multi-variable model should be selected with caution to avoid being affected by overfitting and model complexity. From the VIF value in Table 4, with the increase of the number of variables, the collinearity between variables also gradually increased, which indicates that three is a stable number of variables. Low VIF values and partial correlation coefficients implied low multicollinearity among the independent variables for the population density estimation. In the stepwise process, the higher adjusted R 2 , lower RMSE, and lower MAE (Table 5) were attributed to the addition of POI and NTL data, corresponding to our assumption. Compared with the population density estimation model based on IS data, the method obtains a smaller error [3].
The population density estimation model has a small number of samples: only 141 township samples were used for modelling. Therefore, the model applies to similar scales, such as municipal study areas or small regions. The accuracy of the model was assessed by a 10-fold cross-validation method, indicating that the model was effective in Hefei in 2018. The verification of the model can be further deepened, such as trying to verify it in other urban areas, or applying the model in historical years, and then using census data to verify it. This will be our future research direction. Figure 11 depicts the residuals of the optimal model. We have tested the residuals of these predictions of the multi-variable model. The Moran' I was −0.016, and the p-value was 0.466. Therefore, the spatial distribution of residuals does not conform to spatial autocorrelation. The township with the largest residual is Yafu township, with a residual error of 2602 people/km 2 . Yafu township is located in the southeast of Chaohu ( Figure 1). As can be seen in Figure 3, the ISs are distributed here. Therefore, in future research, we will try to establish population estimation models by calculating the height of urban buildings, which might be very useful to improve the estimation accuracy of the model. Xiaoyaojin Township is located in the southeast of Luyang District, Hefei, covering an area of 2.92 km 2 . The population density here is 23.3 × 10 4 people/km 2 . The northeastern part of Xiaoyaojin Township is dispersed with vegetation, and other parts are full of buildings, as shown in Figure 8a. The north-eastern part is sparsely populated. The building area has high IS and NTL data values, along with a high density of POI data with a concentrated population distribution. Figure 9b cannot distinguish between the north-eastern part and other parts of the country's population distribution differences. The pixels from the northeast of Figure 8c are dark in colour, indicating a sparse population. The others are yellow and orange, indicating a high population density. Therefore, the pixel map ( Figure 8c) can effectively reflect the population distribution in a high-density area.
Located at the junction of the main city of Hefei and Feixi County, Shangpai covers an area of 121 km 2 , with a population density of 2228.29 people/km 2 . From the high-resolution image (Figure 9a), the central area of Shangpai Township is located in the northern part, and several dwelling units are sprinkled in the south-eastern area. Other areas are mainly forests or croplands. The population density of the Shangpai town, displayed in Figure 9b, does not consider the internal population distribution differences. The yellow pixels shown in Figure 9c correspond to the building area in Figure 9a, which represents the area where the population is concentrated.
Xiage Township belongs to Chaohu City (a county-level city in Hefei), covering an area of 184 km 2 , with a population density of 306.27 people/km 2 . It is located on the northern bank of Chaohu Lake, with beautiful scenery. The township (Figure 10a) mainly consists of mountainous terrain with vast tracts of forests, and the residential areas are distributed in the southern part, as shown in Figure 10c.

Discussion
Testing on extracting the proportion of ISs by LSMA obtained a satisfying result, corresponding with the previous findings [29,61]. The map of the ISs showed that their distribution was related to the population (Figures 3 and 7). However, estimating the population with a sole variable may led to a coarse result. Thus, POI and NTL data were considered as added independent variables. Note that the number of independent variables in the multi-variable model should be selected with caution to avoid being affected by overfitting and model complexity. From the VIF value in Table 4, with the increase of the number of variables, the collinearity between variables also gradually increased, which indicates that three is a stable number of variables. Low VIF values and partial correlation coefficients implied low multicollinearity among the independent variables for the population density estimation. In the stepwise process, the higher adjusted R 2 , lower RMSE, and lower MAE (Table 5) were attributed to the addition of POI and NTL data, corresponding to our assumption. Compared with the population density estimation model based on IS data, the method obtains a smaller error [3].
The population density estimation model has a small number of samples: only 141 township samples were used for modelling. Therefore, the model applies to similar scales, such as municipal study areas or small regions. The accuracy of the model was assessed by a 10-fold cross-validation method, indicating that the model was effective in Hefei in 2018. The verification of the model can be further deepened, such as trying to verify it in other urban areas, or applying the model in historical years, and then using census data to verify it. This will be our future research direction. Figure 11 depicts the residuals of the optimal model. We have tested the residuals of these predictions of the multi-variable model. The Moran' I was −0.016, and the p-value was 0.466. Therefore, the spatial distribution of residuals does not conform to spatial autocorrelation. The township with the largest residual is Yafu township, with a residual error of 2602 people/km 2 . Yafu township is located in the southeast of Chaohu ( Figure 1). As can be seen in Figure 3, the ISs are distributed here. Therefore, in future research, we will try to establish population estimation models by calculating the height of urban buildings, which might be very useful to improve the estimation accuracy of the model. In our future research, we will try to combine all independent variables to estimate the population. Additionally, data sources are also important. Geographical and economic data are popular in population estimation studies [62,63], but are difficult to combine with the population because of the spatial gap. Moreover, the lowest scale of economic data is the county (or district) scale, obstructing a deep study on the township scale or smaller. In contrast, the remote sensing data in our study can correspond with the township scale and have obvious advantages over these data types, which are not only easy to access but are updated periodically.
The optimal model developed in this study-incorporating IS data with POI and NTL data in modelling population density and producing a high-resolution map of the population distribution-can effectively and rapidly estimate populations, and has considerable potential in the era of big data. Our future research will try to estimate populations from a smaller scale, such as using sentinel satellites to extract impervious surfaces, in order to obtain a higher resolution of population spatialization results.

Conclusions
Based on IS data, this study gradually introduces NTL and POI data and establishes an optimization model for population density estimation. Exploring the correlation and multicollinearity between variables, the results meet the requirements of establishing a Figure 11. Residuals (on a log scale) between the estimated population density and the corresponding actual population density at the township scale.
In our future research, we will try to combine all independent variables to estimate the population. Additionally, data sources are also important. Geographical and economic data are popular in population estimation studies [62,63], but are difficult to combine with the population because of the spatial gap. Moreover, the lowest scale of economic data is the county (or district) scale, obstructing a deep study on the township scale or smaller. In contrast, the remote sensing data in our study can correspond with the township scale and have obvious advantages over these data types, which are not only easy to access but are updated periodically.
The optimal model developed in this study-incorporating IS data with POI and NTL data in modelling population density and producing a high-resolution map of the population distribution-can effectively and rapidly estimate populations, and has considerable potential in the era of big data. Our future research will try to estimate populations from a smaller scale, such as using sentinel satellites to extract impervious surfaces, in order to obtain a higher resolution of population spatialization results.

Conclusions
Based on IS data, this study gradually introduces NTL and POI data and establishes an optimization model for population density estimation. Exploring the correlation and multicollinearity between variables, the results meet the requirements of establishing a population density estimation model. The 10-fold cross validation of the four models shows that the multi-variable model achieves the best prediction effect. The parameters of the optimal model are applied to 500 m pixels, and the population density distribution map at the pixel scale is obtained after correction. Overall, the multi-variable (including IS, NTL, and POI data) model can effectively predict the population density in interval years and for areas lacking census data.

Data Availability Statement:
The data presented in this study are available upon request from the corresponding author.