Population Spatialization by Considering Pixel‐Level Attribute Grading and Spatial Association

本文研究区域为湖北省武汉市，其街道级行政区划及人口密度等级如图 1所示。武汉市下辖13个区，185个街道，总面积8 569.15 km²。2015年武汉市户籍人口达829.26万人，常住人口达1 060.77万人。13个下辖区中包含7个主城区，即洪山区、青山区、武昌区、汉阳区、硚口区、江汉区和江岸区，占武汉市总人口的61.67%；6个远城区分别为新洲区、江夏区、蔡甸区、黄陂区、东西湖区和汉南区。武汉市不仅具有人口分布众多的主城区，也包含地理范围广阔、人口密度较小的远城区，人口分布情况非常复杂，因此，选择武汉市作为研究区域对于人口空间化研究具有借鉴意义。

Figure 1.  High, Middle and Low Population Density Streets in Wuhan City

本文使用NPP/VIIRS（NPOESS preparatory project/visible infrared imaging radiometer）夜光数据、高德POI及武汉市地理国情普查建筑区数据作为研究数据，详情如表 1所示。夜间灯光数据能反映人类活动，是人口空间化建模的理想数据源^[24-25]。POI数据具有语义丰富且与人口分布高度相关的优点，常被用于人口建模^[14-16]。地理国情普查建筑区数据提供的高精度建筑区轮廓及类别有助于修正夜间灯光溢出的影响，区分不同用地类型，进而辅助人口估算。由于武汉市2015年的POI数据缺失，本文选取2017年数据代替，其他数据源采集时间均为2015年。

PAG-SA的计算与验证流程如图 2所示，共由4个部分组成，包括数据预处理、特征提取、模型训练与估算及格网人口分配。数据预处理阶段对多源数据进行坐标转换、栅格数据重采样、格网信息统计及街道信息汇总。特征提取阶段融合建筑区轮廓数据、夜光数据及POI数据，生成训练及估算时的特征向量。模型训练与估算阶段使用随机森林模型，输入构建的特征向量，输出格网人口权重。格网人口分配阶段基于无房屋无人口原则^[26]约束格网人口权重，在区县级别进行权重归一化并按权重分配格网人口。

Figure 2.  Workflow of Calculation and Validation Process for PAG-SA

数据预处理主要包括多源数据空间基准统一及信息统计。使用ArcGIS等软件将前述数据进行坐标转换，基于不同的格网尺度分别统计格网信息，具体包括：（1）对夜光数据进行坐标转换及重采样，统计各个格网的夜光值；（2）对建筑区轮廓数据进行坐标转换，基于JTS（java topology suite）统计各个格网的建筑区面积占比；（3）对POI进行坐标转换；（4）计算街道人口密度。

PAG-SA综合使用3种特征提取方法构建训练及估算向量：（1）以建筑区轮廓为约束统计格网内各类POI数量特征，以建模不同类型建筑区与POI组合方式对人口密度的影响；（2）使用格网属性分级方法提取基于建筑区类别约束的夜间灯光分级特征，利用格网属性分级方法减小模型跨尺度误差，结合建筑区类别约束缓解夜光值溢出问题；（3）统计邻近格网POI在当前格网中心的核密度估计值，从而建模邻域POI与人口密度间的关系。

作为一种易获取的地理空间数据，POI具有语义丰富且与人口分布高度相关的特点^[14]。基于POI数据进行人口空间化，相比于土地利用类型数据能够更好地保留人口空间分布的细节信息。目前POI数据在人口空间化中应用广泛^[14-16]，但大多只考虑了格网内POI的绝对数量，忽略了建筑区对POI的潜在空间位置约束。不同建筑类型具有不同的人口密度，分布于不同建筑区类型内的POI对人口的吸引力也存在差异。如分布于高密度多层建筑区内的POI对人口的吸引力可能比位于低密度低矮建筑区内的POI更高。因此，PAG-SA基于建筑区类别及空间位置约束，统计不同建筑区类别内的POI数量，以便模型拟合其与人口密度间的相关关系。

具体步骤如下：（1）将建筑区数据和POI数据进行空间叠置分析，统计各个格网中分布于各类建筑区类别内的各类POI数量。假设建筑区种类数为$\mathrm{T}\mathrm{y}\mathrm{p}{\mathrm{e}}_{\mathrm{b}\mathrm{u}\mathrm{i}\mathrm{l}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{g}}$，POI种类数为$\mathrm{T}\mathrm{y}\mathrm{p}{\mathrm{e}}_{\mathrm{P}\mathrm{O}\mathrm{I}}$，则空间叠置后的特征维数为$\mathrm{T}\mathrm{y}\mathrm{p}{\mathrm{e}}_{\mathrm{b}\mathrm{u}\mathrm{i}\mathrm{l}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{g}}\times \mathrm{T}\mathrm{y}\mathrm{p}{\mathrm{e}}_{\mathrm{P}\mathrm{O}\mathrm{I}}$。（2）统计街道内所有格网中各类特征的平均值作为模型训练的输入。

夜间灯光数据能反映居民点、交通道路及产业结构等多种信息，但存在夜光值溢出的现象^[24]，从而影响人口空间化的精度。针对上述问题，有学者提出使用土地利用数据进行约束^[24]，统计格网内城镇用地和农村居民用地的总面积，如果总面积大于0则表示该夜光值有效。该方法能够缓解夜光值溢出问题，但无法应对灯光来源的复杂性。为此，PAG-SA在格网尺度对建筑区数据和夜光数据进行属性分级，将数量信息转换成类别信息，使用类别合并的方法实现数据融合。具体步骤如下：

1）格网属性分级。本文采用自然断点法对夜光值和建筑区类型的面积占比进行分级，通过戴维森堡丁系数（Davies-Bouldin index，DBI）确定最佳分级数量，流程如图 3所示。

Figure 3.  Flowchart of Attribute Grading Based on Natural Breaks

自然断点法是一种数据分级算法，算法原理是对分类间隔加以识别，实现类间方差最大、类内方差最小。DBI系数是一种评估聚类算法优劣的指标，取值范围为[0，+∞），DBI系数越小，表明等级内距离越小，等级间距离越大，计算方法如下：

式中，$\mathrm{D}{\mathrm{B}}^k$表示自然断点法分类数为$k$时对应的DBI系数；${\overline{\alpha }}_x\mathrm{和}{\overline{\alpha }}_y$分别是分级结果中第$x$和第$y$个等级的类内平均距离；${\sigma }_x$和${\sigma }_y$分别是第$x$和第$y$两个等级中心间的距离。

2）特征向量构建。对于一个格网单元，根据分级结果确定各类数据的等级。若格网属于第$t$类数据第$k$个等级，就在其对应的特征向量编码处将特征值赋为1，在第$t$类数据其他等级编码处赋为0。例如，若经过格网分级，夜光亮度值介于[0，100）、[100，200）和[200, 255]分别为第1、2和3类，则对应特征向量编码分别为[1, 0, 0]、[0, 1, 0]和[0, 0, 1]。根据该方法获取特征向量后，将夜光数据与建筑区数据的特征向量按照与运算进行融合以构建组合向量，如下所示：

式中，$\mathit{y}$为组合向量；$\left[l_{\mathrm{1}},l_{\mathrm{2}}\dots l_m\right]$和$\left[b_{\mathrm{1}},b_{\mathrm{2}}\dots b_n\right]$分别表示该格网夜光和建筑区数据的特征向量；$\&$表示与运算。例如，当$l_{\mathrm{2}}$和$b_{\mathrm{3}}$取值均为1，则$l_{\mathrm{2}}\&b_{\mathrm{3}}$等于1，表示该格网夜光亮度为第2等级且建筑密集程度为第3等级。

在街道单元尺度，若第$i$个街道的总格网数为$N_i$，属于第$t$类建模数据第$k$个等级的格网数量为$N_i^{t,k}$，则其特征向量可表示为：

前面提取的特征只包含格网本身的POI语义信息，导致格网间的空间关联被人为割裂。为此，本文使用核密度估计提取邻域POI特征，以刻画邻近格网POI对当前格网人口分布的影响及距离衰减效应^[27]。核密度估计是分析点事件分布和识别热点^[28-29]的一种常用方法。相关研究及本文实验表明，核函数的选择对结果影响不大^[14]，带宽（搜索半径）是核密度估计的主要参数^[22]。

本文选用密度函数较为平滑且使用场景广泛的高斯核作为核函数：

点对象核密度估计方法^[29]如下：

式中，$k$为核函数；$\mu \left(s\right)$表示位置$s$处的核密度估计值；$r$为带宽；$d_{is}\mathrm{为}i$点到当前位置$s$的距离。本文针对每一类POI，通过比较多种带宽取值获取相对最优带宽。

本文对POI进行核密度估计时，位置$s$为当前格网中心点，如图 4所示。图 4中红色点表示当前格网中心，蓝色点表示POI，圆柱半径表示搜索半径。各POI相对于当前格网中心的核密度估计值记作$k(d_{is}/r)$，格网中心最终核密度估计值为$\mu \left(s\right)$。统计各个街道内所有格网的平均值作为模型训练的输入。

Figure 4.  Illustration of Kernel Density Estimation upon Grid Cell Centre

本文选用随机森林构建回归模型，其原因为：（1）该模型对异常值和噪声具有较高的容忍度^[13]。在人口空间化中，由于数据源的多样性及人口分布的复杂性，在特征向量中往往存在异常值却难以发现。如以街道数据进行训练时，面积较小的街道可能存在数据分布极端的训练样本。（2）融合多源数据导致特征维数增多，造成筛选及降维困难，而随机森林模型能够处理高维数据，避免人为特征选择。（3）随机森林模型中决策树相互独立，利于并行实现，训练速度快。

PAG-SA的训练与估算流程如图 5所示，训练阶段输入街道级别特征拟合街道人口密度，估算阶段输入格网级别特征生成格网人口权重。

Figure 5.  Workflow of Training and Estimation for PAG-SA

针对人口权重可能存在无建筑区但权重非零的问题，本文根据无房屋无人口原则^[26]，将无建筑区的格网人口权重赋0。

经过建筑区数据修正后，对各个区县进行权重归一化并将区县人口按照权重分配到各个格网中，根据格网与街道的映射关系计算街道人口。格网人口计算方法为：

式中，$i$表示第$i$个区县；$j$表示第$j$个格网；$\mathrm{S}{\mathrm{I}}_{ij}$表示第$i$个区县第$j$个格网的最终人口值；$\mathrm{S}{\mathrm{I}}_i$表示第$i$个区县的人口总值；$W_{ij}$、$W_{iu}$分别表示第$i$个区县第$j$个和第$u$个格网权重值；$N_i$表示第$i$个区县的格网总数。

本文选取平均绝对值误差（mean absolute error，MAE）、均方根误差（root mean square error，RMSE）、决定系数（R²）3种指标进行精度评价，计算公式如下：

式中，MAE反映人口估算误差的绝对值；RMSE刻画人口估算值与真实值之间的偏差程度；R²度量人口估算值与真实人口的拟合程度；$N$表示街道总数；$\mathrm{P}\mathrm{r}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{c}{\mathrm{t}}_i$表示街道$i$的估算值；$\mathrm{R}\mathrm{e}\mathrm{a}{\mathrm{l}}_i$表示街道$i$的真实值；$\overline{\mathrm{R}\mathrm{e}\mathrm{a}\mathrm{l}}$表示所有街道真实人口的平均值。

为了验证特征提取方法融合多源数据的有效性，将仅使用POI密度的方法、综合使用POI、夜间灯光及建筑区数据进行特征向量直接拼接的方法与PAG-SA在200 m、500 m和1 000 m共3种格网尺度下进行精度对比，结果如图 6所示。图 6中，特征向量直接拼接的方法不考虑特征之间的关联，POI统计值为POI在街道或格网内的密度，夜光统计值为街道或格网内的平均夜光亮度，建筑区统计值为建筑区在街道或格网内的面积占比。

Figure 6.  Accuracy Comparison of Three Feature Extraction Methods

图 6表明PAG-SA能够有效提升精度，且不同格网尺度的效果存在显著差异。直接拼接方法在3种格网尺度下均出现R²下降、MAE/RMSE上升的情况，说明使用建筑区及夜光数据在街道尺度上训练构建的模型对格网尺度不适用，引发模型跨尺度问题。同时，不恰当的数据融合方式可能导致精度降低。而相对于仅使用POI的方法，PAG-SA的拟合优度及准确度在3个尺度下均有所提升。MAE下降4%~16%，且200 m格网尺度精度最优，随着格网尺度增大，精度逐渐降低，说明不同尺度下特征的精细程度和表达能力不同。

为了进一步验证特征提取各步骤的有效性，本实验对比了4种特征提取方法，实验结果如图 7所示。方法1仅使用POI密度作为参照；方法2提取建筑区类别及空间位置约束的POI特征；方法3在方法2基础上提取夜间灯光分级特征；方法4在方法3的基础上引入POI核密度特征。

Figure 7.  Accuracy Improvement at Each Feature Fusion Step of PAG-SA

由图 7可知，除了1 000 m格网外，各特征提取步骤均有助于精度提升，但是不同格网尺度的效果不同。（1）从格网尺度上看，200 m的R²呈现接近线性的上升趋势，而500 m和1 000 m的R²先上升后趋于平缓，说明PAG-SA相对适合较小格网尺度下的数据融合；（2）使用建筑区类型及空间位置约束的格网POI特征在3个格网尺度下均取得显著的精度提升，原因是该方法有助于提升POI特征的精细度；（3）融合基于建筑区类别约束的夜间灯光分级特征在200 m下有一定的精度提升，而其他两个尺度下提升较小，说明格网尺度较小时，夜光分级特征能更真实地反映人口分布规律；（4）邻域格网POI特征提取在200 m尺度下精度提升较明显，而500 m提升较小，在1 000 m甚至出现精度下降。原因在于尺度较大时格网自身已包含相对丰富的信息，引入邻域格网特征反而增大误差。

PAG-SA在200 m格网尺度下的人口空间化结果如图 8所示，其人口分级采用自然断点法。从空间分布模式上看，武汉市人口呈现中心城区集聚且周边多核的空间结构，人口值大于275的格网主要分布于中心城区。

Figure 8.  Results of PAG-SA in 200 m Grid Size

对于不同人口密度区域，其人口建模特征的空间分布存在差异，空间化结果亦呈现不同模式。为此，本文使用自然断点法将武汉市185个街道按人口密度值划分为高、中、低3个密度等级区域，并将PAG-SA与WorldPop、GPW及中国公里网格人口数据集进行对比。街道人口密度分级结果如图 1所示，其中低密度区60个街道，中密度区69个街道，高密度区56个街道。各街道误差如图 9所示，其绝对误差为人口估算值与人口普查值之差，相对误差为绝对误差与对应街道人口普查值的比值，高低估街道定义为相对误差绝对值大于0.1的街道。图 9中估算人口散点图中蓝、绿、红3种颜色分别表示低密度区、中密度区及高密度区。

Figure 9.  Comparison of Accuracy Between PAG-SA and WorldPop, GPW and PopulationGrid_China

由图 9可知：（1）从3种评价指标的数值上看，PAG-SA在高、中、低3种密度区域相对于对比人口数据集均有更小的误差。（2）从拟合效果上看，PAG-SA的散点大体集中分布于对角线两侧，而对比数据集的散点较分散且距离对角线较远，说明PAG-SA能够更好地拟合真实人口分布。（3）从绝对误差空间分布上看，总体而言，GPW与中国公里网格人口数据集在高、中、低3种密度区域均存在大量绝对误差大于40 000的街道，WorldPop的绝对误差主要分布于中密度区域。而PAG-SA除武汉市东南区域及其他零星分布街道外，均有较为明显的精度提升，显著降低了误差等级。（4）从相对误差上看，PAG-SA的高估和低估街道数量相对均衡，而WorldPop与中国公里网格人口数据集在低密度及中密度区域易高估，而在高密度区域易低估。原因是WorldPop及中国公里网格人口数据集的估算范围较广，提取的特征不够精细，因此估算结果较为平均化。GPW数据集在高、中、低3个密度区域均易高估，原因是GPW数据集基于格网内的行政单元面积进行人口分配，未顾及其他组合因素的影响。（5）PAG-SA在武汉东南区域没有明显的精度提升，原因是该区域POI等设施较齐全但实际人口较少，使用武汉全部街道训练构建的模型不适用于该区域，可通过分区域训练提高精度。

考虑到各类POI的辐射范围^[22]不同，为了获取核密度相对最优带宽，实验对比13种带宽下（200~1 000 m间距为100 m，2 000~5 000 m间距为1 000 m）单独使用每一类POI构建邻域格网特征时的精度。由于带宽选择受到格网大小的影响，本文实验对比了200 m、500 m及1 000 m这3种尺度下带宽选择对精度的影响。实验表明，200 m格网下3种精度评价指标对带宽选择最敏感，且200 m下所获取的相对最优带宽区间包含500 m及1 000 m对应的最优区间，因此本文采用200 m尺度格网开展带宽选取实验。带宽选取实验结果如图 10所示，图 10中虚线和实线分别表示引入POI核密度特征前后对应的精度评价指标。

Figure 10.  Impact of Bandwidth on Estimation Accuracy for Different POI Types

由图 10可知，不同类型POI的相对最优带宽不同，且各类POI在各自相对最优带宽处均能提高精度。（1）与引入核密度前对比，在某些带宽区间使用核密度后MAE/RMSE降低，同时R²提高，说明合适带宽下各类POI的核密度特征对提高精度均有效。（2）从MAE/RMSE的变化趋势上看，科研教育、住宿及金融服务这3类POI的整体变化幅度不大，对带宽的选择不敏感，而医院、休闲娱乐、购物及餐饮4类POI形成了明显的波峰。因此，选择核密度带宽时需结合POI类型，不同POI类型设置同一个带宽可能会引入误差。（3）从相对最优带宽上看，医院和休闲娱乐的相对最优带宽约为3 km，可能原因是医院和休闲娱乐场所的辐射距离较大，空间服务范围较广；而餐饮、购物及居民小区的相对最优带宽较小，说明这3类POI总体上辐射距离较小，空间服务范围较为有限。科研教育、住宿及金融服务的带宽取值对精度的影响不敏感，反映出这几类设施的辐射能力较强且空间服务范围广泛的特点。各类POI的相对最优带宽取值范围如表 2所示。

本文提出一种顾及格网属性分级与空间关联的人口空间化方法：（1）基于自然断点法分别对建筑区密度和夜间灯光值分级并融合二者构建组合向量，在行政尺度使用各等级格网占比信息作为训练输入，以减小模型跨尺度误差；（2）通过实验为不同类型POI选择合适的核密度估计带宽构造邻域格网POI特征；（3）基于叠置分析统计建筑区类别及空间位置约束的格网POI特征，以便刻画多种属性不同空间聚合方式与人口密度间的关联关系，提高特征建模精细度。以武汉市为实验区域，通过与WorldPop、GPW及中国公里网格人口数据集的对比验证了方法的有效性。实验结果表明，本文方法街道尺度MAE远小于对比数据集，在高、中、低人口密度区域均取得较好的拟合优度，并有效提升空间化精度。同时本文特征提取方法中各步骤的有效性均得到验证，且200 m格网尺度精度提升最为明显。不同类型POI的辐射作用范围不同，合理的核密度带宽阈值与各类POI的社会职能相关，需通过实验选取。

本文方法存在以下不足有待进一步研究。本文利用建筑区轮廓和POI间的空间关系进行数据融合，提取的特征虽然保留POI在不同建筑区类型内的数量信息，但未考虑空间分布信息。POI的空间分布模式，如均匀、随机或聚集，可能对人口分布产生影响进而有助于刻画人口分布^[30]，今后将尝试引入空间分布特征。同时，本文使用枚举方式选择POI核密度带宽，今后可研究最优带宽的自适应提取方法^[31]以提升带宽选择的效率及可解释性。由于依赖于建筑区轮廓及POI等细粒度数据，因此本文模型无法直接迁移到相关数据缺失的区域。但其特征建模方法依然具有一定适用性及参考价值，如格网属性分级、核密度估计及空间叠置分析等，在后续工作中尝试将上述建模方法迁移到不同类型的区域进行验证与分析。