CN111080501A - 基于手机信令数据的真实人群密度时空分布估算方法 - Google Patents

Abstract

本申请提供了一种基于手机信令数据的真实人群密度时空分布估算方法。人群密度时空分布估算的方法包括：使用深度卷积生成对抗网络(DCGAN)训练数据，生成虚拟人口密度分布；获取虚拟实验中的手机信令数据的呼叫详细记录(CDR)；绘制估算的人口密度分布。本申请的估算方法不仅可以使政府对公共资源管理变得容易，而且可以使在手机上计算显示这种的实时人口密度分布成为可能。

Description

基于手机信令数据的真实人群密度时空分布估算方法

技术领域

本发明涉及基于手机信令数据的真实人群密度时空分布估算方法。

背景技术

人口动态分布对城市和交通规划，公共安全预警，灾害影响评估[7-9]和流行病学建模存在着深远的影响。然而，在研究人的活动时，以更高(的时空分辨率超过12小时)估算和绘制更详细的人口分布仍然是一个挑战。具有较高时空分辨率的人口分布可用于许多应用中，例如，人口分布预测可以为地方当局提供基础数据和参考，以避免或减轻踩踏事件；在城市尤其是在中国的城市中，有很多共享的自行车为市民服务，政府可以通过根据人员分布动态地配置资源来提高其利用率。另一种商业用途情况是在食品外卖市场中，如果提供者知道人口分布在何时何地以更高或更低的密度分布，他们可以动态而有效地安排商品和劳动力资源以获取更多利润。

最近，许多研究人员和公司越来越多的使用手机数据的进行实际应用，并且通过挖掘大数据的呼叫详细记录(CDR)来估计和映射人口动态分布。当我们需要动态的人口分布来应用于其他科学研究时，CDR作为一种基于当前100％手机渗透率的位置数据，成为了的最为流行的数据集。而GPS数据由于涉及到用户隐私的问题一般都难以获得，尤其是在大范围的空间中。为了满足以更高的时间分辨率(例如每小时或每分钟)进行人口动态分布的要求，我们可以考虑基于CDR数据进行估算。

但是，由于CDR中的手机用户不能等于实际人口，因此在将其转换为真实人口密度分布的过程中需要解决许多问题。例如，人们通常在白天比晚上打电话的频率更高，因此生成并存储在CDRs数据集中的记录更少，因此，如果我们仅使用移动电话用户或记录号来代表实际的人，晚上不使用手机的用户将不会被记录，之后就会导致直接用CDR估算人群分布的时候产生极大的误差。因此，现有一些研究侧重于如何将手机数据转换为真实的人口密度分布的估算方法，除了静态分布和耗时的自身缺点外，没有真实数据基础作为对比也成为其估算精度存在不确定性的重要原因。

发明内容

为了可以测试估算方法的精度，本申请提出了一个虚拟实验，在其中生成虚拟人口密度分布并进行估算工作，最终以通过对比测试本申请的估算结果精度。

本申请提供了一种人群密度时空分布估算的方法，包括：使用深度卷积生成对抗网络(DCGAN)训练数据，生成虚拟人口密度分布；获取虚拟实验中的手机信令数据的呼叫详细记录(CDR)；绘制估算的人口密度分布。

在上述方法中，其中，所述呼叫详细记录包括用户与基站的交互时间、基站的位置以及用户的加密ID。

在上述方法中，还包括：通过建立虚拟基站分布并且基于所述虚拟基站分布来获取所述呼叫详细记录。

在上述方法中，其中，人口密度分布的公式如下所示：

其中，Na为区域的人数，Ns为从用户的加密ID中随机抽取的用户的数量，M为Na除以Ns的倍数，Ni为区域中的每个Voronoi多边形内的总和，i代表每个Voronoi多边形的ID。

本申请的有益效果如下：

1)本申请提出了一种简单高效的方法，可以在其他方法无法返回实际动态结果的情况下，根据CDR数据有效地估算动态实际人员的密度分布；

2)使用由深度学习方法DCGAN生成的虚拟人群密度分布的方法建立，首次以通过建立虚拟实验的方法来测试本申请的估算方法的鲁棒性和准确性；

3)本文提出的估算方法可以提供快速高效的绘图过程，以小时时间分辨率绘制实际人口密度分布图，可用于了解人们的动态热点分布和政府公共资源的管理；

4)此外，随着移动GIS(地理信息系统)不断的发展，本申请的算法使这种可在手机上计算的实时真实人口密度分布成为可能。

附图说明

图1示出了本申请的方法的流程示意图。

图2示出了本申请的基于腾讯定位大数据的DCGAN训练图像举例：(a)至(c)分别来自北京，南京和上海的城市，时间为下午20:00，日期为2019年4月29日。

图3示出了基站移植采样区域及其Voronoi多边形分布。

图4中的(a)和(b)分别示出了整个北京和基站采样区域中Voronoi多边形的面积分布结构。

图5示出了通过DCGAN生成人口密度分布的灰度结果图。

图6示出了鉴别器(d_loss)和生成器(g_loss)的损失(loss)变化图。

图7示出了在一个示例测试中，通过不同的采样点具有不同的渲染结果的估计人群密度分布的密度图。第一个是原始基础地图，第二个是相同的分布，但是密度已分为15类，显示了图像和KDE渲染的栅格网格图之间的差异。当本申请从1000到10,000个移动电话用户的不同点进行采样时，其他10张图像是根据本申请的方法估算的结果。

图8示出了相似度与不同采样数的输入图像之间的关系。

图9示出了输入图像分辨率与最后的相似度之间的关系。

图10示出了基站用户交互频率在一天中的变化。

图11示出了北京市2015年2月17日的高精度时空人口密度分布。

图12示出了顺义区和首都国际机场区位图。

图13示出了北京市顺义区2015年2月17日的估计人口密度分布：该图第一列中的4张图像分别是上午1：00、1：30、2：00和2:30的分布，而第二列中的下4张图像是从09:00到10:30AM，最后一张是当天最后4个半小时的图像。

图14示出了2010年人口普查北京市人口密度分布图。

图15示出了三种分布图结果的精度对比：(a)显示三个比较的RMSE：估算和普查数据，用户和普查数据，记录和普查数据；(b)显示与普查数据相比的估计结果；(c)说明了单用户使用结果与单记录使用结果的比较。

具体实施方式

下面的实施例可以使本领域技术人员更全面地理解本发明，但不以任何方式限制本发明。

在大空间中获取人口密度的传统方法是人口普查，该方法虽然较为准确，但由于其周期(通常每10年)过长，因为只能用于人口密度的静态分布。另外，在空间分辨率上，人口普查仅限于县或镇等区域，无法再精确到小区或建筑物。作为研究人口分布的一种手段，人们已经积极探索使用遥感技术。但是，仅使用遥感数据只能在粗略/中等分辨率下得出映射的人口密度，这通常限制了此类数据在小规模尺度上的适用性。之后，将遥感技术和其他数据源结合起来进行估算的方法受到了许多研究的关注。另外，本申请可以使用视频监视或图像来估计获得较小规模区域的人口。但是，它仅适用于较小的空间，不适用于城市规模。因此，当GPS数据由于隐私原因而无法访问时，使用手机大数据在更详细的时间和较大的空间范围内动态估计和映射人口是一种趋势。

由于在如此大的空间尺度上没有更高的时间分辨率的动态实际人口分布的地面真实数据，因此有必要使用可以产生模拟人口分布的方法建立虚拟实验来验证估计方法的准确性和稳定性。在小规模的室内房间和住宿社区中，基于代理的建模(ABM)是一个很好的系统，可用于人群踩踏的疏散和交通流量管理。另一项研究概念化并指定了一种称为LUDAS(土地使用动态仿真器)的多主体系统，该系统代表农村森林边缘的人地景观耦合系统，以探索替代方案以改善生计并减轻土地利用的负面影响变化，从而支持土地利用规划中各个利益相关者之间的谈判过程。除土地利用外，气候变化也被视为影响人口流动的关键因素，并且该文说明了如何结合基于计划行为理论的基于主体的模型来预测基于证据的未来移民变化以响应未来人口，经济，社会和气候变化。

近些年来，生成对抗网络(GAN)框架不断的发展。它主要是原理是通过对抗过程来估计生成模型。其中有两个模型：生成模型G，用来捕获数据分布，以及区分模型D，用以估计样本来自训练的概率数据，而不是G已经过训练。这种深度学习模型可以训练图像类型并生成伪图像。之后，有学者引入了具有一定体系结构约束的一类称为深度卷积生成对抗网络(DCGAN)的CNN，以弥补CNN成功实现监督学习和基于GAN的无监督学习之间的差距。GAN和DCGAN经过测试，对生成伪人脸图像，查看图像甚至音乐具有良好的效果。作为具有地理特征的特殊图像，本申请可以在DCGAN中训练实际人口密度分布图，来生成虚拟的人口密度分布。

数据

本申请的内容分为两个主要部分，其工作流程如下(图1)。第一部分是本申请的实验内容，旨在生成虚拟人口密度分布，然后展示本申请提出的估算方法如何工作；第二部分是本申请的方法在北京市中的应用。在第一部分中，本申请有四个步骤：第一步是使用DCGAN训练腾讯公司的定位数据所绘制的人口密度分布，之后并生成虚拟人口密度分布。第二步，本申请在实验中随机抽取一定的样本点作为虚拟实验中的CDR。然后，本申请的方法用于绘制估算的人口密度分布；在最后一步中，本申请使用余弦相似度方法比较这两个分布。在第二部分中，本申请将所提出的估计方法应用于中国北京的CDR数据，绘制高精细度时空人群密度分布并做了进一步分析。

腾讯公司的微信一直是中国最受欢迎的社交软件，根据《2018年微信年度数据报告》有超过10.1亿用户登录使用微信。根据《凯度中国社交媒体影响力报告2018》，在北京，上海和广州等中国较大城市中，微信的普及率至少达到了96％的公民。除微信软件外，腾讯公司还开发了许多其他移动应用程序，例如社交应用程序QQ，腾讯视频和QQ浏览器，它们也以较高的普及率而流行。所有这些应用程序都需要其用户打开定位权限，以便一旦本申请安装了其中之一，腾讯公司便可以通过其智能手机GPS获取用户的位置。根据用户位置数据，腾讯提供了一个实时平台来显示其用户的热图，并且密度分布可以通过GIS软件轻松绘制例如ArcGIS。根据腾讯应用程序的普及率和用户数量，腾讯定位大数据将更接近中国，特别是大城市中实际人口密度分布的真实数据。

尽管腾讯定位大数据更接近真实数据，但由于其只适合渗透率高的中国城市，因此并不能作为替代人口普查的估算方法。但是作为可以训练生成虚拟人群分布的数据集，本申请可以将其利用在本文中。本文首先基于2019年4月29日至30日在北京，南京和上海的每小时定位腾讯大数据的分布，提取了随机固定范围20km*20km空间范围内的6000个分布图。本申请显示了渲染之后具有相等密度值的15级灰度图像后的分布(样本如图2所示)。这些数据集用于训练模型，然后通过DCGAN生成模拟的人口密度分布，然后进行虚拟实验。

在方法部分的申请中，本申请使用了来自中国移动运营商公司(北京)的匿名个人呼叫详细记录(CDR)，其中包含有关国际移动用户识别码(IMSI)的信息。IMSI是每个SIM卡或HLR的国际唯一代码/VLR识别网络上的用户；timestamp是记录交互式事件何时发生的时间戳；CI记录时的小区标识对应的基站位置。该数据集是匿名收集用于科学研究的，每当用户呼叫或发送短信时，都会记录一条新的数据记录。本套数据集包括2015年2月17日(星期四)在中国北京的超过3亿用户的48亿条记录。

表1.CDR结构

本申请使用了全球移动通信系统(GSM)工程参数内部结构(表2)中所示的所有51216个移动基站的位置。但是，由于基站数量非常多，以至于所确定的纬度和经度是相同的，因此本申请将重复定位的基站从51216到17445进行了组合。每个移动基站的覆盖范围可以近似为Voronoi多边形，当用户使用电话进行呼叫或发送短信时，可以通过在电话所连接的特定移动基站的范围内找到其位置。

表2.全球移动通信系统工程参数

根据特征CI，IMSI和时间戳，每个记录中的用户CI将由IMSI识别，并将相应的基站位置视为用户的位置。本申请称CDR中时间戳的准确性为记录的时间分辨率，而在此数据集中，时间分辨率为1秒，并且CDR每1小时(3600秒)存储在一个逗号分隔值(CSV)文件中。在空间上，由于城市中手机的使用量较大，因此基站的密度也较大，一个站覆盖半径为100-500m的圆形区域。而在郊区，受需求和地形的影响，单站覆盖较大的圆形区域，半径约为500m至1000m。

方法

第1步：基于DCGAN构建虚拟实验环境

为了测试本申请的分布估算方法，本申请需要将其对比基础分布。但是，由于没有真实的基础分布数据，本申请建立了一个虚拟实验，在其中生成虚拟人口的分布，然后使用本申请的方法进行对于虚拟CDR的估算与绘图。在此部分中，本申请使用腾讯定位数据集作为训练数据然后构建DCGAN来为下一个流程生成虚拟分布。因为本申请的方法没有考虑不同时间的影响来估计人口分布，所以在比较过程中不需要考虑人群活动随时间的变化。因此，在这一步骤中本申请仅考虑具有人口分布格局特征的静态分布。

训练分布图来自北京，南京和上海的腾讯定位数据来自，如图2所示。由于人口的分布在空间上具有相似的分形特征，因此本申请将中国的这三个典型城市进行混合训练，以生成用于虚拟实验的人群分布。然后，本申请遵循Radford,A.；Metz,L.；Chintala,S.Unsupervised representation learning with deep convolutional generativeadversarial networks.arXiv preprint arXiv:1511.06434 2015中DCGAN的原始集，设置学习率为0.0002，批量大小为64×64，训练数据集是提取的6000张地图，输入图像大小设置为108×108像素(灰度)。本申请训练了具有1000个时刻点的采样地图，最终生成了64个虚拟分布图像。在下一步中，本申请从64个结果中选择了10个，以对本申请的估算方法进行交叉检验。

对于每个选定的实验图，将根据图像上的灰度颜色生成随机点。例如，如果本申请定义该区域中有50万人，则根据每个颜色区域成比例地将相同数量的点放置在地图上。本申请将最终生成的分布图定义为基础分布。

第2步：随机采样构建CDR

本申请将已记录在CDR中的手机用户作为实际人群的随机样本，因此，在此步骤中，将提取基线分布中特定人数的点作为手机用户。

本申请还需要建立虚拟基站分布，这是生成虚拟CDR的基础。为了使虚拟实验更接近真实情况，本申请选择从北京提取的具有基站分布的实际区域移植到测试区域，以及实验中人口数量和实际人口数量的比例(本申请将其设置为Pexperiment)应该等于测试区域与整个北京区域的比例。同时，基站数量也需要近似等于比例。

在本申请中，该区域的面积是一个边长为20公里的正方形，它是通过在空间中随机选择搜索的，直到包含425个基站(425/17445等于400/1441万平方公里，1441万是整个北京的区域)如图4所示。所有Voronoi多边形的区域分布结构如图4所示，可以说明两个区域的相似性，从而使实验更接近实际情况。现在，基准虚拟人口分布和基站分布已构建完成。

第3步：高精细度人口密度估算方法

本部分将具体介绍本申请的估算方法的细节。在介绍之前，本申请需要报告的该方法所基于的事实条件。由于某些城市的人口数量多年来增长或减少缓慢，而本申请需要的是每小时甚至每分钟的细粒度时空密度分布，因此一年中的变化可以忽略。尽管在大城市中，尤其是那些可以吸引大量游客的城市，但随着常住人口的增加，几天内的人口变化相对稳定。例如，根据中国国家移民局在2015年第一季度进行的移民和边境检查的综合统计，从2015年1月1日到2015年3月31日，从北京出境的人数为62,321,614，前往北京时为62,447,662。这两个数字确实很接近，根据北京市统计局发布的结果，2015年北京市常住人口为2170.5万人，北京每日总人口的平均变化率仅为0.5％。因此，如果只考虑人口数量而不考虑诸如年龄和性别之类的个人信息，则人们的活动可以被视为在封闭区域内发生的过程。在这种情况下，尽管它们始终一直从一个位置移动到另一位置，但仅在封闭区域内。

基于上述城市规模上人口稳定的前提，本申请开始介绍这种估算方法的细节。首先，基于研究区域是具有固定人数Na的封闭区域的前提，本申请在假设进行研究的每个时间点上，根据CDR从IMSI代码中随机抽取相等数量的Ns个移动用户。之后，计算每个基站多边形中的手机用户数量(如图1所示)。第三，本申请使用该区域的总人数Na除以Ns得到常数M，然后将从CDR中提取的每个手机用户的数据集乘以M，同时将每个基站的M个点随机放置在地图上区域。最后，将核密度估计(KDE)方法用于估计总体密度分布的输出。

在本申请的研究中，本申请将记录呼叫或发短信事件的基站位置视为移动电话用户的位置。Thiessen多边形算法用于为每个移动电话基站创建Voronoi多边形，以便定义一个区域内电话用户的位置(图1)。根据上一段中的分析，基于城市规模中人口稳定的前提，本申请认为手机用户与基站的交互是随机事件，这意味着手机用户只是拥有CDR记录的随机抽样人员。另外，由于智能手机的高普及率，采样过程比较平均。因此，本申请可以使用样本人员活动来估计整个人员活动，并将移动电话用户的分布视为缩小的实际人员分布比例。就时间尺度而言，如果本申请想获取某个时间段内的估计人口分布，则可以提取该时段内具有空间身份特征(CI)的手机用户，以代表实际人群的样本数据。

为了运行KDE，地理信息系统软件(例如ArcGIS)计算了适当的搜索半径，以得出KDE的估计值。在不同的空间分辨率(例如800米)下的密度，以生成由100m*100m像元的正方形像元组成的栅格网格。

在这里，KDE用于执行自动搜索以探索事件分布的热点，并使用复杂的距离衰减来测量事件密度的变化。基于GIS的KDE估计方法主要使用移动窗口来计算和输出每个网格单元的点或线密度。给定样本(x1，x2，…，xn)，是一个独立的，均匀分布的样本，该样本是使用Rosenblatt-Parzen核估计计算得出的，在点x，f(x)上具有分布密度函数f的簇：

其中k()是内核函数；h>0是内核的平滑参数，称为带宽；(x-xi)是从估计点x到采样点xi的距离。

在处理KDE时，带宽h的确定或选择对计算结果有很大影响。随着h的增加，空间中点密度的变化更加平滑，但是密度的结构被掩盖了。当h减小时，估计的点密度变化在Voronoi多边形之间可能会突然变化。在特定应用中，有必要根据物理环境中不同土地用途(例如，堤岸或火车站的存在)测试不同的h值，以便探索可以匹配实际情况的核密度表面。

KDE估计的具体步骤是：1)定义搜索半径，以计算落入带有滑动圆的圆内的事件的数量；2)根据密度精度要求确定输出栅格大小；3)通过核函数计算每个事件对循环域中每个网格的密度贡献；4)将每个栅格的密度值分配给每个事件在栅格搜索半径中的密度贡献值；5)输出每个栅格网格的密度值。在这项研究中，本申请使用ArcGIS10.5内核密度分析工具确定北京地区800m以内的人口密度分布的特征作为默认半径，空间分辨率为每单元100米。

第4步：虚拟实验中实际值与估算值的比较

每次使用KDE并为其网格中的每个像元设置密度值后，都将输出栅格网格。经过上述过程之后，本申请将在实验中以栅格栅格格式获得估算的人口密度分布。由于基线分布是由DCGAN预测的图像生成的，因此，为了比较基线分布和估计的分布，需要将栅格网格转换为具有相同标准的图像，这意味着在符号体系中渲染器的颜色或灰度分配系统ArcGIS中的过程(可从https://pro.arcgis.com/en/pro-app/help/mapping/layer-properties/symbolize-feature-layers.htm获取，于2019年1月3日检索)应与网格。生成的基线图像以108×108像素的灰度输出，以使图像中的单元格可以轻松表示地理单位，除了本申请以200×200的分辨率将估计分布的网格转换为灰度外，还将尺寸调整为相同的大小，其中每个像素代表一个边长为100地理单位米的正方形区域，而灰度值表示每平方公里单位人口的密度。

然后本申请使用余弦相似度方法来测量两个图像的相似度值，公式如下：

其中S是两个图像A和B的相似度；Ai表示图像A的第i个值的值，而Bi表示图像B的第i个值的值，这是本申请中的浓度值。将图像A和B中的像素值转换为矢量后，通过计算累积点积与这些2范数的比率，可获得最终的余弦相似度结果。S的最大值为1，在该范围内，该值越小，图像之间的区别越大。

结果与讨论

虚拟实验

第1步：生成虚拟人口密度分布

从DCGAN方法生成结果的64个结果(如图5所示)中，本申请随机选择其中10个进行10次虚拟实验。本申请比较了10对图像的相似性，以显示本申请方法的性能。最后，本申请将该方法应用于整个北京，并绘制了一天中每半小时的动态估计人口密度分布图。

在DCGAN过程中，鉴别器损失量化了鉴别器能够区分真实图像和伪造图像的能力。它将辨别器对真实图像的预测与一个1s数组进行比较，并将辨别器对虚假(生成的)图像的预测与一个0s数组进行比较。生成器的损失量化了其欺骗鉴别器的能力。直观地，如果生成器运行良好，则鉴别器会将伪图像分类为真实图像。图6显示了当本申请使用DCGAN训练并生成模拟分布时，鉴别器损耗和发电机损耗随历时增加而变化的细节。生成器损耗在开始时波动，然后在历时达到约300时达到约0.3的稳定点，这意味着鉴别器倾向于将伪造图像检测为真实图像。鉴别器的损失从减少开始，并且在历时达到300时也保持在约1.7的稳定性。这种损耗收敛通常表示该模型达到了最佳状态，在此过程中它改进得很多，已经学到了足够的知识。实验的良好效果显示所生成的是可靠的虚拟人群分布图。

第2-3步：估算结果

为了对图像数据的统一对比，本申请将所有图像渲染为灰度图：白色区域表示较高的值，黑色区域表示较低的值。得到10个基础人口密度分布后，下一步是根据密度分布将随机点放入地图中。在本文的上一节中，本申请将500,000个点放入相应的地图中，并使每个测试的密度分布与图像相同。一个点意味着一个人，本申请随机抽取样本点作为与基站发生交互的手机用户，并且该事件已记录在模拟CDR中。

为了检测从虚拟人员个数的采样数与估算结果之间的关系，本申请将采样数设置为1000到25,000，并绘制栅格密度图像，然后得到10组图像。图7显示了基础分布和相应的估计结果的示例，其中使用了1000至10,000个不同的样本数(本申请总共采样到25,000个)。根据先前的研究，城市中的人们聚集在一起呈放射状和环形分布，例如中国上海。当这些特征显示在本申请的结果中时，也已经检验了该结论。最后，本申请将每幅图像与每10次测试中的基本分布的灰度图进行比较。

第4步：对比结果

图8显示了10组平均相似度随采样数增加而变化的规律。当输入图像所对应的人员点采样数从1000到12,000时，平均相似性指数急剧增长，而当采样数大于12,000后，其值稳定在0.982左右。本申请还计算了当抽样数量为500,000(与实验人数相等)时的相似性指标，这意味着本申请将所有生成的数据用于本申请的方法进行比较。10组的平均比较结果为0.985，这意味着基于余弦相似度方法，本申请的方法的最高准确性为98.5％。但是，在实际情况下，不可能获得100％的样本，因为在上一节中已说明的原因(例如，并非所有移动电话用户都在特定时间段内与基站进行交互)。综上所述，当从12,000个开始增加采样数时，效率是非常接近的，占整个测试数据的41.7％。

使用余弦相似度方法计算相似度指标时，输入图像的分辨率可以更改。在初始图像中，最大输入分辨率为200*200，而最小输入分辨率为1。通过将输入分辨率从1更改为200(正方形边长的像元数)，可以在相同的条件下计算相似度。

图9说明了不同样本数量1000、5000、10000、15000和20000的相似度和输入分辨率之间的关系。显然，随着输入分辨率的提高，所有相似度都从近似值1逐渐减小，这是由两个图形的所有像素的平均灰度级相同引起的。当分辨率提高到50*50时，相似度指数以不同百分比达到最低点，该百分比由样本数决定，但所有行从50*50一直保持稳定直到分辨率200*200。因为在使用KDE时本申请将每个像素的边长设置为100米，所以根据图9，本申请的估算方法的空间分辨率可以高达500米。当每个基站小区的地理长度小于400米时，相似性指标是稳定的；而当长度大于400米时，相似性指标会随着输入分辨率的降低而增加到1。但是，输入分辨率1*1没有意义，尽管输入地理分辨率高于400*400米，但相似度不会增加。这意味着当使用KDE输出最终密度图时，像元大小可以设置为400米，而这浪费了机器内存和使用范围，导致栅格边的长度比地理400米短。

不同样本数的估计结果说明了该方法的不同水平能力。同时，输入图像的分辨率会影响准确性。在下一部分中，本申请将在北京的情况下使用具有特定样本数量和KDE输出地理单元大小的方法，并使用手机数据绘制某一特定日期半小时的人口密度分布图。

应用实例：北京市高精度时空人口密度分布估算

CDR中手机用户的提取与分析

作为中国的首都，北京在2015年2月17日的CDR数据集将作为本申请方法的应用。根据上文的分析，本申请可知在2015年第一季度，即2015年1月1日至2015年3月31日，从北京出门的人数为62,321,614，而去北京的人数为62,447,662。这两个数字很接近，根据北京市统计局发布的结果，2015年北京市常住人口为2170.5万人，北京总日均人口的平均变化率仅为0.5％。因此，本申请认为北京是一个封闭的城市，每天的人口数量稳定，而没有考虑每个人的特征。

应用本申请的方法时，应首先从CDR中提取样本人口，此过程不考虑发生频率而提取唯一移动电话用户。但是，移动电话用户使用电话的频率不是平均的，因为在夜晚，例如在2:00AM，使用电话的人数明显减少，而在白天，例如10AM，则有更多的人使用电话。电话使用他们的电话进行呼叫或发短信，该设备与记录事件的基站进行交互，因此本申请将事件称为站与用户的交互事件。然后，计算2月17日每1小时前30秒钟的CDR中的手机用户，然后本申请通过公式f＝N/t计算频率，其中N是手机唯一用户数，t是30秒。基站-用户交互频率与一天中时间之间的关系如图10所示。很明显，从一天开始到7:00AM，较少的人打来电话或发短信，而在7:00AM则急剧增加。以及上午10:00。接下来的8小时会有一些波动，并且从下午18:00开始迅速下降直到这一天结束。

条形图(图10)说明，不同的采样持续时间和不同的时间会导致不同的精度。例如，在凌晨2:00的第一秒中，CDR中记录的人口只有70,000，这意味着如果本申请要估计这一秒的实际人口密度分布(在这种情况下，采样持续时间为一秒)，公式3.1中的M为310(北京总人口数21,700,000除以70,000)。记录的用户数占总数的0.3％，这对应于图9中的采样数1.5，并且响应相似度相对较低。但是，如果本申请从凌晨2:00开始在一个小时内提取出全部唯一的移动电话用户，到凌晨3:00，移动电话用户的总数约为1000,000(由于连续两秒之间存在大量重复用户，因此它远远少于3600乘以70,000的总和)，这对应于图9中的24,000个样本现在通过本申请的方法进行估算的准确性更高。相反，在上午10:00当天，在10:00A.M的第一秒的CDR中记录了24万人，公式3.1中的M为90(北京总人口21,700,000除以240,000)。记录的用户数占总数的1.1％，对应于图9中的5.5个采样，相似度相对较低。同时，如果本申请从上午10:00起的一个小时内提取出全部唯一的手机用户，到凌晨11:00，移动电话用户的总数比一秒钟的数字要多得多，因此估算的准确性更高。因此，为了在更多的手机用户样本中达到更高的准确性，采样持续时间从凌晨2:00开始。

在本申请的研究中，整天的手机用户抽样数量固定为100,000，这意味着抽样率即北京的手机数量与总人口之比为0.4％，与图8中对应的采样数为2000时的0.94相似度对应。基于采样率，最高采样分辨率可能是凌晨2:00的2分钟30秒在上午10:00，但是，为了以规则的时间间隔绘制动态估计的人口密度分布，每半小时的前2分钟提取手机用户进行绘图。尽管本申请可以以更精细的时间分辨率(例如一整天2分钟)绘制，或者以更高的估计精度绘制那些采样更多手机用户(但精细度较低的时间分辨率)的采样，但本申请只是展示了本申请方法的效率通过快速估算过程以及一天中空间尺度的明显变化。

半小时时间人口密度分布的动态估计

在2015年2月17日提取北京手机用户样本后，利用相应的CI一天中每半小时识别用户的位置48次。按照第3.2节中的过程，由于该部分的样本数为100,000，而整个北京的人口数为2,170,000，因此，在属于该多边形的区域内随机生成了代表每个Voronoi多边形中的人口的点数的217倍。然后用KDE方法对点分布进行了48次映射，绘制出了像元大小为400×400平方米的密度分布栅格网格。最后，按照凌晨00:00的第一个网格。在2015年2月17日的最大和最小密度值上，本申请通过几何间隔方法(可从https：//pro获得)将所有值分为32个分类，并用蓝色至红色的颜色条渲染了网格，红色表示密度较高，而蓝色表示密度较低。其他47个网格使用与第一个网格相同的符号系统方法，所有这些都构成了估计人口密度分布图像数据集。

图11显示了北京在2015年2月17日全天以每2小时较高的空间分辨率估算的人口密度分布。显然，在一天中的24小时内，密度分布每次都在不同区域动态变化。在空间规模分布上，在北京中部的东城区，西城区，人的密度始终高于密云区，怀柔区等郊区。在每个行政区域中，该区域中心镇的人口密度分布都比其他区域居民高得多。例如，在北京北部怀柔区中心的怀柔镇，人们的密度要高得多，而在周边大多数城镇，如怀柔镇以北的雁栖镇，人口密度要低得多，原因之一是这个小镇上有很多人居住的山区。就时间尺度上的分布变化而言，人们在夜间(例如2:00AM)更分散，而在白天(例如11:00AM)则更紧密地聚集，这是因为大多数工作场所都分布在然而，在北京市中心，该地区的住宅用地较少，许多人不得不以较低的房价在郊区购买公寓。因此，夜间人口密度较高，而白天则较低，这是郊区人口分布的特征。相反，夜间的密度较低，而白天的密度较高，这是城市地区的特征，因为工作日的主要工作场所是工作日。

进一步鲁棒性检测：将结果与CDR中的单个用户/记录数进行比较

尽管没有地面事实数据可用来验证估计结果是否准确，但可以将每个行政区划的人口普查数据用于比较夜间(00:00至8:00AM和22:00至12点)的人口密度分布：00PM)，因为这段时期的人是相对静态的。在本申请中，本申请使用2010年的北京人口普查数据来比较CDR的分布，图14显示了北京人口密度2010年分配给322个行政单位。

本申请从2015年2月17日上午00:00每2小时开始提取1分钟CDR。同样的过程也应用在上午8:00，下午22:00上，最终形成了6个数据集。对于每个数据集，本申请将每条记录视为一个人来绘制第一个分布，提取记录在CDR上的所有手机用户来绘制第二个分布，最后使用本申请的方法来估计和绘制第三个分布。然后，根据北京市的行政单位，本申请从每个单位内的三个分布图计算平均密度。下一步是将这三个组分别与人口普查分布进行比较。图15中的(b)所示的小提琴图的左侧部分显示了322个管理单位的估算结果的平均密度值，红色部分是普查数据。比较这两个密度值的一部分，很明显，它们非常相似，并且在所有人口普查数据都相同的6个时间点对估计的一个密度值具有稳定性。根据人口普查数据，在322个单位中的平均人口密度为每平方公里8402.954人，而在00：00、02：00、04：00的估计人口密度分别为7307.027、684.985、6640.536、7032.84、7763.540和7368.336，2015年2月17日06：00、08：00和22:00。但是，由于本申请提取的CDR较少，因此从单个记录编号和具有位置特征的用户编号绘制的密度远低于人口普查。尽管它表明从短时间内提取的CDR数量不足以使用这两种传统的绘制方法，但是如果目标是在短时间内以较少的计算内存绘制人群分布，则本申请的估计方法是有效的。

同时，图15中的(a)展示了三个分布结果的RMSE结果和六个时间点的人口普查数据，还显示了夜间人们活动的准确性和真实性，其中RMSE约为5000用于估计方法，对于其他两个分布，大约为13500。

在包括城市和农村地区在内的大空间尺度上对实际人口分布动态估计和绘制具有深刻而实际的应用。除了用于城市和交通规划之外，细粒度的动态人口分布还可以用于公共安全预警，灾难影响评估和流行病学建模，从而使政府，商人和公民受益。为了满足以较低的价格和较短的时间在细粒度但大范围的时空尺度上获得人类活动的需求，并在将来将其整合到手机中，本申请提出了一种估算方法，该方法使用以下方法绘制动态人口密度分布图手机数据。

为了评估本申请方法的准确性，本申请的估算方法已在由深度卷积生成对抗网络(DCGAN)生成的虚拟情况下的10个实验中使用。DCGAN的结果表明，使用该网络基于地面真实数据集或类似的数据集(如腾讯定位大数据)生成人口密度分布的效果非常好，这意味着在以后的研究中，研究人员还可以使用DCGAN生成虚拟分布图。估算结果表明，随着参数的变化，余弦相似度指标有所不同，但总体上具有较高的估计能力。最后，作为本申请方法的应用，以北京为例，以2015年2月17日为例，对半小时动态人口密度分布进行了估计和绘图，以分别分析了整个北京和和顺义区的人口时空分布特征。

本申请提出的估算方法适用于在研究时间内人群流量稳定的大城市，并且适用于对时间和空间上分辨率高的条件情形下。每个网格单元的最高空间分辨率为400*400平方米，而时间分辨率可以在几分钟到几小时之间变化，这取决于估计精度的不同要求。较高的估计精度要求需要更多的移动电话用户样本数量，这意味着较低的时间分辨率，而较低的估计精度要求则需要较少的移动电话用户采样数量，这意味着较高的时间分辨率。总而言之，更高的精度和更高的时间分辨率无法同时使用。但是，较高的时间分辨率(例如导致较低精度的分钟数)将足以绘制人口密度分布的动态变化，分析活动热点并检测人群分布特征。

最后，由于人口普查数据可以更加真实表现夜间人口分布情况，本申请又基于简单的将用户作为真实人口，基于简单的将记录作为真实人口进行了直接密度绘制，然后分别和本申请的方法结果进行了比较。结果还显示了本申请的估计方法的准确性，真实性和内存节省，该方法的RMSE约为5000，而其他两个分布超过13500，该过程节省了资源，内存和时间来提取和分析CDR数据，但精确度还很高，这就是本申请将本申请的方法视为快速而准确的方法的原因。

未来进一步的工作将在虚拟实验部分的步骤2中考虑不同空间中CDR中手机用户的差异。这种考虑将使随机采样的人点随着CDR的复杂化而误差降低。例如，在一个工作日晚上有许多白领加班的商业区域，人们的活动将比住宿区域丰富，这可能导致手机用户的CDR中的记录分布不均，但在本申请的实验中是随机抽样的没有考虑。此外，由于本申请的方法需要较少的内存和对中央处理单元的计算请求，因此移动设备(例如手机)可以将该方法作为功能强大的移动GIS软件集成到移动应用程序中，这可能是有用的功能，并且对移动GIS的发展有着实际意义。

本领域技术人员应理解，以上实施例仅是示例性实施例，在不背离本申请的精神和范围的情况下，可以进行多种变化、替换以及改变。

Claims (4)

1.一种人群密度时空分布估算的方法，包括：

使用深度卷积生成对抗网络(DCGAN)训练数据，生成虚拟人口密度分布；

获取虚拟实验中的手机信令数据的呼叫详细记录(CDR)；

绘制估算的人口密度分布。

2.根据权利要求1所述的方法，其中，所述呼叫详细记录包括用户与基站的交互时间、基站的位置以及用户的加密ID。

3.根据权利要求1所述的方法，还包括：通过建立虚拟基站分布并且基于所述虚拟基站分布来获取所述呼叫详细记录。

4.根据权利要求2所述的方法，其中，人口密度分布的公式如下所示：

其中，Na为区域的人数，Ns为从用户的加密ID中随机抽取的用户的数量，M为Na除以Ns的倍数，Ni为区域中的每个Voronoi多边形内的总和，i代表每个Voronoi多边形的ID。

Images