本发明属于城市规划设计领域，特别涉及一种基于手机信令的城市中心评估方法。
背景技术：
从整体上评估城市中心体系是城市空间结构研究的重点之一，也是城市总体规划的基础工作之一。过去的研究一般采用普查统计数据，例如人口普查数据、经济普查数据。由于普查数据基于一定的普查单元，空间精度较差，识别的仅仅是可能包含中心的街道等统计单元。因此，城市中心体系的识别与评估工作受到很大制约。
技术实现要素：
本发明的目的是提供一种基于手机信令的城市中心评估方法。
一种基于手机信令的城市中心评估方法，包括预处理、识别和评估这三个步骤，
预处理又包括以下步骤：
a1，根据手机信令给出的位置信息，选择工作日每天0:00-6:00作为识别居住地的时间段，在该时段内的所有位置点中选择与其它位置点距离的和最小的点作为当天居住地，再以同样的方法识别白天(10：00-16:00)的稳定就业地，最后得出具有稳定居住和工作地的手机用户和手机用户数量；
a2，将每个手机用户每天的空间位置数据按照临近关系进行归并，距离500m内归并为一组，认为是其一次出行目的地；
a3，统计该手机用户去过每个目的地的频率和时长，将去过天数小于等于两天，且停留时长大于2小时的出行作为生活性出行，提取其中双休日出行od(起点—终点)记录。“o”来源于英文origin，指出行的出发地点，“d”来源于英文destination，指出行的目的地。od的含义是交通起止点。
根据od记录，识别两类城市居民行为，一是通勤，二是生活性出行，通勤目的地集中的地区，即潜在的就业中心；生活性出行目的地集中的地区，即潜在的生活中心，通勤和生活性出行od信息作为识别与评估两大城市中心的基础数据；
a4，基于识别的od信息，利用现有技术中的空间插值算法获得密度曲面，将识别的通勤就业目的地进行插值，获得就业岗位空间分布密度曲面，用来识别就业中心，
将识别的生活性出行目的地进行插值，获得生活性出行目的地的人次密度曲面，用来识别生活性中心，密度曲面均采用gis的栅格数据形式存储，
识别步骤包括：
将城市中心的定义为：要素分布密度显著高于周边区域、且对周边要素分布具有一定影响的空间单元，
具体计算形式是采用简单线性回归，对每一个目标栅格单元i周边r范围内的单元j，计算单元j到单元i的距离dij，并提取j处的访问密度值yj，对yj、dij进行一元线性回归：
yj＝ω1×dij+ω0+ε，j∈r
对合适半径r的选择，借鉴贪心算法，事先认定必须识别为中心的点，然后从较小的半径逐步扩大搜索半径r，每一次搜索半径都会产生一次识别结果，从这些结果中选择事先认定的中心识别效果最好的，作为最终识别结果。实验表明，r取值在1200m左右效果较好。
拟合获得ω1、ω0、决定系数r2三个系数，以ω1<0且r2>0.1作为识别中心的边界的依据，以0.1作为决定系数门槛，考虑相关系数在0.3以上认为具有中等相关性，
评估步骤包括：
c1，对识别出的中心，对平均到访人次密度进行k均值聚类分析，来划分中心的等级，将k值定为5，对应“市级中心”、“市级副中心、”“地区级中心”、“重点镇”、“一般镇”规划实践中的5个等级；
c2，通过汽车交通时耗评估生活中心的可达性，行车速度大致采用道路设计行车速度乘以0.8折减系数，用来减少与实际车速的误差，分别计算各栅格单元到达高等级生活中心、全部生活中心所需要的时间；
c3，通过势力圈逐级分析考察生活中心服务的人口规模，势力圈即中心地主要服务的空间范围，在该范围内居民前往对应生活中心获得服务的概率最大。
本发明利用手机信令数据识别并评估城市中心体系的方法，基于由手机信令数据识别的出行起讫点信息识别城市的就业中心体系、生活中心体系，并通过可达性分析、势力圈的方法评估。手机信令数据是手机用户通话、短信、切换基站、位置更新等行为的被动记录，可通过数量庞大的基站位置连续追踪手机用户时空轨迹，具有大样本、全时性等特点，且空间精度远高于一般统计数据。
附图说明
图1是本发明的方法流程框图。
具体实施方式
如图1所示，本发明的基于手机信令的城市中心评估方法，包括预处理、识别和评估这三个步骤。
对于手机信令数据预处理，以上海市连续两周共14天匿名手机信令数据为例。数据处理借鉴已有方法，包括三步：首先，选择工作日每天0:00-6:00与其他记录点的距离和最小的点作为当天居住地，再对工作日(10天)的居住地重复距离和最小法识别稳定居住地(对于多天识别居住地差异过大的用户，认为其没有稳定居住地)，最后得出具有稳定居住和工作地的用户752万人，占同年上海常住人口的31％。第二步，将每个用户每天的空间位置数据(居住地1000m之外)按照临近关系进行归并，距离500m内归并为一组，认为是其一次出行目的地。第三步，统计该用户去过每个目的地的频率和时长，将去过天数小于等于两天，且停留时长大于2小时的出行作为生活性出行。该阈值中，“2小时”是因为手机周期性更新信令的周期约2小时；“2天”是根据14天内用户在居住地和工作地之外的任一地点重复出现的天数分布情况，选择与生活性出行较为贴近的“2天及以下”部分。最后，提取其中双休日出行总计约443.9万条od记录作为识别与评估中心的基础数据。
预处理过程识别了两大类城市居民行为，一是通勤，二是生活性出行。通勤目的地集中的地区，即潜在的就业中心；生活性出行目的地集中的地区，即潜在的生活中心。经过预处理的通勤和生活性出行od将作为一下识别与评估两大中心体系的基础数据。
基于识别的od信息，利用已有空间插值技术获得密度曲面。将识别的通勤就业目的地进行插值，获得就业岗位空间分布密度曲面，用来识别就业中心。将识别的生活性出行目的地进行插值，获得生活性出行目的地的人次密度曲面，用来识别生活性中心。密度曲面均采用gis的栅格数据形式存储。
传统的城市中心识别技术同样基于普查数据。由于普查数据基于一定的普查单元，例如乡镇街道单元，因此空间位置精度低，无法反映城市地块尺度的活动分布。在此制约下，传统中心识别技术一般包括两个步骤：
(1)先将普查单元视为点，采用局部加权回归(locallyweightedregression)来获得平滑的要素密度曲面，例如就业岗位的空间分布密度栅格，此时需要选择平滑系数(指参加局部回归的观测值的个数与所有观测值个数的比例)；
(2)然后筛选出那些密度值远大于拟合的密度平面上对应的平滑值的区域，作为候选的城市中心。常用做法是选择残差值在5％显著水平显著为正的空间单元作为备选的城市次中心；
(3)最后，还需要检验这些候选中心对周边点密度值的影响的显著性。即，将周边点的密度值作为因变量，周边点到市中心的距离、周边点到候选点的距离作为自变量，构建一个复杂的检验模型，检验到候选中心的距离是否对周边点的密度值有影响。最后通过拟合显著性来判断候选中心是否作为最终识别中心。
在这一系列过程中，规划师需要面临：选择平滑系数、选择置信区间、选择检验模型形式、选择检验模型中全局密度分布函数、甚至还要添加一些虚拟变量来辅助计算。这些选择稍有不同就会对结果产生很大影响，因此传统方法的可复制性差，在不同地区规划项目中应用的可行性极低。另外，由于传统普查单元数量少(上海只有5000多个普查区)，而手机基站数量极大(本发明实施例采用的有3万余个)，如果将传统传统方法应用到手机信令数据上会使计算量大大增加，远超一般计算机运算能力。
对于识别，本发明将城市中心的定义：要素分布密度显著高于周边区域、且对周边要素分布具有一定影响的空间单元。具体计算形式是简单线性回归，对每一个目标栅格单元i周边r范围内的单元j，计算其到单元i的距离dij，并提取j处的访问密度值yj，对yj、dij进行一元线性回归：
yj＝ω1×dij+ω0+εj∈r
对合适半径r的选择，借鉴贪心算法，事先认定必须识别为中心的点，然后从较小的半径逐步扩大搜索半径r，最后选择认定的中心识别效果最好的结果。结果表明r在1200米左右效果较好。拟合获得ω1、ω0、r2(决定系数)三个系数。以ω1<0(“要素分布密度显著高于周边区域”)且r2>0.1(“对周边要素分布具有一定影响的地区”)作为识别中心的边界的依据。以0.1作为决定系数门槛，考虑了相关系数(决定系数的平方根)一般在0.3以上认为具有中等相关性。本发明首先确定一些公认的商业中心，然后在多个阈值下对比这些商业中心的识别效果，发现在0.1-0.2之间差异不大，且0.1效果更好。
该方法本质上是在识别局部曲面与理想的“山峰状”中心的相似程度，任何一个栅格单元只要它周边曲面满足该条件，相当于也满足了中心的理论定义，即可纳入中心的范畴。需要指出的是，该方法仅适用于空间精度较高、连续性较好的数据，例如本文采用的手机信令数据。该方法的优点是：编程过程简单、计算速度快、适合大范围大数据量的城市区域环境。共识别了185个生活中心和191个就业中心。
传统的普查数据如人口普查、经济普查数据，受统计单元精度所限，经插值后的密度曲面无法反应客观的分布状态，不一定适用本发明的方法。
对于评估，包括步骤：
(1)对识别出的中心，本发明提出对平均到访人次密度进行k均值聚类分析，来划分中心的等级。将k值定为5，对应“市级中心”、“市级副中心、”“地区级中心”、“重点镇”、“一般镇”等规划实践中的5个等级。
(2)借助arcgis已有工具，通过小汽车交通时耗评估生活中心的可达性，行车速度大致采用道路设计行车速度乘以0.8折减系数，用来减少与实际车速的误差。分别计算各栅格单元到达高等级生活中心、全部生活中心所需要的时间。
(3)本发明提出通过势力圈分析法考察生活中心服务的人口规模。势力圈即中心地主要服务的空间范围，在该范围内居民前往对应生活中心获得服务的概率最大。由于势力圈的分析过程需要在中心之间进行比较，而等级差别太大的中心之间比较就失去了意义，因此采取逐级比较分析。
对于本发明的城市中心识别效果检验，包括：
(1)总体检验。识别中心185个，总面积仅占市域3.6％，但集聚了30.5％的生活出行量。识别就业中心191个，总面积仅占市域4.3％，但集聚了28.6％的就业岗位。显然，识别结果显示出了城市中心具有的高度集聚的基本特征，总体上效果较好。
(2)选择识别出的中心里形状较为不规则的同时进行检验，与城市实际功能分布进行对比，验证是否识别结果存在问题。一共选择了6个生活中心，包括“四平路控江路”、“延吉中路”、“共康”、“闵行体育公园”，以及两个与规划边界出入较大的中心，即“陆家嘴张杨路”、“中环(真北)”。
“四平路控江路”因主要公共设施呈l形分布而形成了心形中心；“延吉中路”、“共康”则是由于主要公共设施沿主干道、地铁线呈线性分布形成了长条形中心；“闵行体育公园”则由公园、市场形成两个中心组合式的杠铃形态。“陆家嘴张杨路”规划的商业中心边界沿世纪大道呈带状，识别的中心位于八佰伴和世纪大道地铁站之间，集聚了大量商业广场、综合体、酒店、文化中心等功能，东西两端则以商务办公楼宇为主，显然识别结果准确反映了实际情况；“中环(真北)”规划的商业中心位于真北路立交旁的百联购物广场地块，但是在西侧地块也有着数个商业综合体、商业广场以及大型公共服务设施(沪西医院)，识别结果这两个中心基本呈对称形态，准确反映了实际情况。
(3)高等级生活中心的功能检验。
采用微博poi对识别的生活中心进行功能检验。数据采集自2015年12月新浪微博网站，微博poi具有功能信息较为准确、空间覆盖较为全面的特征，可反映某个地区的功能及其比重、检验生活中心识别结果。统计识别的高等级生活中心内部的poi功能信息，如表1所示，可见高等级生活中心都提供了综合的生活服务，与规划经验较为相符，表明识别结果准确有效。
表1高等级生活中心的功能统计
对于评估方法效果检验，包括：
(1)通过可达性分析，统计了各个乡镇街道到达高等级生活中心和所有生活中心的平均时间，然后将时间可视化制图，即可快速获得对生活中心空间分布的评估结果。这样的快速评估是传统数据和传统调研方法所无法企及的。
对于高等级生活中心，大多数街道单元都能在20分钟内到达；仅有宝山北部、浦东东部、闵行东南部以及远郊的较多乡镇到达时耗在20分钟至30分钟；大于30分钟的有浦东曹路镇和远郊边缘地区的乡镇。对于全部生活中心，大多数街道单元都能在10分钟内到达；存在问题的主要是远郊乡镇，尤其是青浦西部和浦东东南部，到达时间在20分钟以上。
(2)通过势力圈分析，快速获得各个等级的中心服务的人口规模，可作为总体规划布局中心的依据。如表2，对生活中心势力圈的分析可见，地区级生活中心服务总人口最多，是城市生活服务的最主要提供者；市级生活副中心的平均服务人口达120万人，分别是地区级生活中心的约2倍、郊区重点镇的约4倍，恰好呈现梯度倍数关系。这些分析结果也是传统数据和方法无法精确获得的，而规划实践迫切需要的。
表2高等级生活中心服务人口规模统计
基于本发明提出的识别方法，制作了基于arcgis10.2平台的中心识别工具。该工具无需安装即可运行于arcgis10以上平台。
工具是用来根据一个要素分布曲面识别中心，对于其它城市大数据也有应用的可行性。当输入一个栅格，工具会在给定的邻域上进行局部一元回归拟合，根据拟合结果判断中心范围。因此，需要用户输入一个栅格，一个起始的邻域半径，工具会在给定的邻域内进行局部拟合。根据拟合结果，通过k、决定系数来识别中心。k<0表示目标单元高于周边邻域，决定系数阈值保证了目标单元显著高于且对邻域要素分布有显著影响力。如果用户输入了检验要素类(即用户根据经验和实际情况给定的，必须识别为中心的点)，工具会从给定的邻域半径开始逐渐减小搜索半径，直到所有给定的检验中心都被识别了出来，分析才会结束。如果有检验中心没有识别出，将在检验中心要素类创建字段予以标注。
工具包括7个输入栏：
“工作空间”：工作空间是用来暂存临时文件、保存最终输出结果的位置。建议选择一个文件夹。请尽量保证路径无中文。
“输入栅格”：用于识别中心的输入栅格。采用较小的像元尺寸可能获得较好的效果。输入栅格的坐标系必须为投影坐标系，单位必须为“米”。
“搜索半径”：用于对目标栅格进行局部回归时的邻域搜索范围。
“检验要素类(可选)”：用于检验中心识别的要素类，是用户确定一定以及肯定属于中心的输入。当输入检验要素类时，建议设置一个较大的搜索半径，如1500-2000米。工具的搜索半径会自动从大(输入的搜索半径)到小(输入搜索半径的一半)逐渐尝试，每一个搜索半径识别出的中心都将与输入的检验要素类对比，若检验的中心全部识别出，则识别过程停止。
“识别中心时的k阈值”：表征目标栅格与周边栅格相对关系的参数。默认k小于0，表示目标栅格大于周边邻域的栅格。
“识别中心时的决定系数阈值”：表征目标栅格对周边栅格值大小的影响作用。
“输出要素类”：输出要素类是一个面类型。输出结果将采用与输入相同的坐标系。