1 引言

城市是空间、功能和人类活动的组合^[1]。人们在城市空间的活动表征了城市空间的功能分异。在时间、位置以及出行目的的共同影响下,城市人群选择性地在城市内部区域之间移动,到达不同的场所,进行不同的活动。海量的交通数据中隐藏着城市人群出行的时空规律,也蕴含着城市功能分区的信息。从海量的交通数据中挖掘城市人群的出行模式,有助于发现城市区域的功能和区域间的相互作用,揭示城市的空间结构。

目前,城市交通数据被用于城市人群出行距离的研究^[2,3,4]、城市人群出行时空规律究^[5,6]、城市热点区域研究^[7,8]、城市空间交互作用^[9,10,11,12]等方面;结合城市POIs数据,交通数据还可用于判断城市功能区的分布^[12,13]。由于交通运营数据量较大,数据隐含的维度也比较复杂,用一维向量或者二维矩阵形式难以完全表达多维度信息。张量模型则可以充分发挥其多维度的特点,从不同的维度描述数据,已经成功应用于信号处理、计算机视觉、数据挖掘等方面^{[14,15,16,17]}。在城市交通出行方面,Sun等^[18]使用张量分解模型探究新加坡城市公交车数据中隐藏的城市人群出行的时空规律,Wang等^[19]通过张量分解模型揭示北京出租车轨迹数据中隐藏的城市区域模式,他们的工作都是将上车点（O）、下车点（D）与其它信息一起构建统一的张量。但是城市居民出行活动通常在一定的距离范围内^[2,3,4],将O点数据和D点数据放在同一个张量中进行模式挖掘,必然使得OD点相互影响,从而导致同一类模式出现地理空间聚集现象。

作为城市流动人群的主要代步工具,出租车的运营状况和上下车点在很大程度上反映出人群移动的时空规律。相对于出租车完整轨迹路线的随机性,出租车OD数据更能反映出乘客的行为特征,并且出租车OD数据包含轨迹移动的空间信息和时间信息,可以反映出租车用户在不同时段中上车点与下车点对应的热点区域,揭示城市在不同时段中热点区域的空间分布与人群集聚现象,表达城市区域的动态性。而目前利用手机通讯、出行轨迹等大数据研究城市功能的研究中大多只考虑静态的区域功能,未反映出区域功能的动态变化^{[12-13,20-21]}。本文利用北京市出租车OD数据,对北京市出租车用户的出行行为与活动进行分析,推测出行目的及其对应的城市空间功能,并结合兴趣点（Point of Interest, POI）的分布,提高推测的准确性。为解决多维度信息的表达和挖掘,本文选择张量模型;为避免相同区域模式的地理聚集现象,本文使用出租车OD数据分别构造上车张量和下车张量,通过张量分解,发现出租车用户上下车行为在不同维度的时空特征,获取不同时段对应的空间模式,实现出租车用户出行的功能区动态识别。

2 研究区和数据源

2.1 研究区域概况

本文实验数据包括北京市六环路网数据、出租车轨迹数据和POI数据,统一坐标系为WGS 1984 UTM Zone 50 N（WGS84椭球体的通用横轴墨卡托投影第50分度带）,研究区域范围：X[422 211.60 m, 475 711.60 m],Y[4 393 462.30 m, 4 447 962.30 m],研究区域在指定投影坐标系下为53.5 km×54.5 km的矩形区域。本实验将研究区域划分为500 m×500 m格网,不仅简化问题模型,而且降低计算的复杂度,划分结果如图1所示,总计11 663个格网,表现为109行×107列,格网编号为从上至下、从左到右依次编号为0-11 662。

显示原图| 下载原图ZIP| 生成PPT

图1   研究区域

Fig.1   The study area

2.2 实验数据与预处理

本研究出租车轨迹数据时间范围为2012年 11月1-16日,从轨迹数据中提取出租车OD数据,删除重复、不完整以及范围超限数据,其中六环区域内有效轨迹为2 531 618条记录,将每一条出租车轨迹的O点与D点映射到格网中,每一条出租车轨迹记录字段见表1。

城市POI数据泛指一切可以抽象为点的地理对象,尤其是与人们生活密切相关的地理实体,可以反映出城市空间区域的功能属性。本文加入北京六环区域的POI数据识别城市功能区特征。本文POI数据通过高德地图开发平台提供的API获取,该平台默认将POI数据分为23大类,263个小类。根据实验目的与需求,本文筛选、拆分与合并不同类别的POI数据,最终获取总计760 200条,总计13个类别（餐饮服务、购物服务、生活服务、体育休闲服务、住宿服务、风景名胜、住宅、政府机构、科教文化、交通设施、金融保险和公司企业）,每条记录包括POI名称、经度、纬度、行政区与类别属性。

3 研究方法

张量是一种多维数组,不同的维度可以具有不同的坐标系统。一阶张量是一个向量,二阶张量是一个矩阵,三阶及以上阶数的张量则是高阶张量,图2是一个三阶张量X∈$\mathbb{R}^{I\times J\times K}$。

显示原图| 下载原图ZIP| 生成PPT

图2   三阶张量

Fig. 2   A third-order tensor

3.1 张量分解

张量分解为奇异值分解和主成分分析的高阶推广,主要包括Tucker分解和CP （Canonical Polyadic Decomposition）分解。其中Tucker分解是一种高阶主成分分析,它把一个张量分解为一个核心张量沿着每一个模乘上一个因子矩阵,每个模上的因子矩阵称为张量在每个模上的基矩阵或者主成分。以3阶张量X∈$\mathbb{R}^{I \times J\times K}$为例,对其进行Tucker分解,分解模数为（P×Q×U）（式（1））,分解结果如图3,其中核心张量G∈$\mathbb{R}^{P\times Q\times U}$,因子矩阵包括A∈$\mathbb{R}^{I\times P}$,B∈$\mathbb{R}^{J\times Q}$,C∈$\mathbb{R}^{K\times U}$。如果P,Q,U分别小于I,J,K,G可以看作是X的压缩,压缩之后的核心张量存储空间明显小于原始张量存储空间,从而达到高维数据的降维,起到数据压缩与特征提取的作用。因子矩阵的每一列为一个主成分,代表该因子特征空间的一种模式,每一行为某一事物在不同模式的值,表示这一事物与不同模式的相关程度,模式值越大相关性越大。核心张量表示因子矩阵之间的相互作用或者联系强度; Xˆ为核心张量与三个因子矩阵的乘积,是一个近似张量^[14]。

X≈X∼=G×A1×B2×C3=∑p=1P∑q=1Q∑r=1Rgp,q,rap×bq×cr=〚G,A,B,C〛（1）

显示原图| 下载原图ZIP| 生成PPT

图3   张量分解-Tucker模型

Fig. 3   Tensor decomposition-Tucker model

3.2 张量分解的计算

张量分解求解问题为：给定一个张量X∈$\mathbb{R}^{I\times J\times K}$,求解核心张量G∈$\mathbb{R}^{P\times Q\times U}$和因子矩阵A∈$\mathbb{R}^{I\times P}$,B∈$\mathbb{R}^{J\times Q}$,C∈$\mathbb{R}^{K\times U}$,使X与 Xˆ差距最小。KL （Kullback-Leible）离散度可以量化相似结果的质量,KL离散度越小,2个相似结果越接近。X与 Xˆ的KL离散度定义如式（2）。

式中：X_{i, j, k}表示张量中的元素值, Xˆi,j,k表示近似张量中的元素值。因此求解问题转化为最优化KL离散度 D(X||Xˆ)问题。

张量分解的核心张量与因子矩阵中的元素 可为正或负,从数学角度观点看,分解结果中存在负值是正确的,但是负值在实际问题中可解释性差,通常是没有意义的,同时张量分解结果不具有唯一性。为解决以上问题,Mørup等^[22]提出SN-TUCKER （稀疏非负Tucker分解）,该模型在张量分解过程中加入非负约束与稀疏约束。对于任意一个非负张量X,SN-TUCKER通过非负约束可以分解出非负核心张量和非负因子矩阵;通过稀疏约束,对核心张量和因子矩阵进行降维,从而减少张量分解结果的歧义性。SN-TUCKER利用乘法更新规则迭代求解,从而保证G,A,B和C的非负性。SN-TUCKER利用1-范数实现稀疏约束,使结果尽可能的稀疏以获取主要特征,从而减少分解结果的歧义性。

4 实验分析

4.1 实验流程

本实验首先进行数据预处理,然后通过张量 分解获取北京六环区域的区域模式。划分到同一种模式中的区域在一天的某个时间段内,具有相似的上下车量,城市居民在相邻的时间段出行或者在相邻的位置上下车,对应的出行模式可能存在一定的相似性,即在某一个时间段内,上车点所在区域或者下车点所在区域具有相似的功能特征;在此理论基础上结合区域活跃时间段和POI分布进行出租车用户出行的功能区特征识别,实验流程如图4所示。

显示原图| 下载原图ZIP| 生成PPT

图4   实验流程图

Fig. 4   The flowchart of experiment

4.2 张量构造与参数设置

4.2.1 张量构造

本实验构造上车张量X_O = （上车点,时段,日）,维数为11 663×24×16,如图5所示,总计16张天切面,每张切面共有11 663行24列,其中每一张垂直切面表示一天;垂直切面的每一行表示一个区域,每一列表示24 h中的一个时段,其中0时段对应 0-1点,其它时段以此类推;垂直切面中格网元素值表示某个区域在某个时段中的出租车上车量。同理构造下车张量X_D = （下车点,时段,日）。

显示原图| 下载原图ZIP| 生成PPT

图5   张量构造

Fig. 5   Tensor construction

4.2.2 张量分解参数设置

本文对张量X进行分解,求解核心张量G,因子矩阵：A区域模式矩阵、B时段模式矩阵和C日模式矩阵。在进行张量分解之前,需要确定张量分解的分解模数和稀疏约束参数。“时段”模数参考其它文献^[4,18-19]和经验确定为4,即早高峰时段、日间时段、晚高峰时段以及夜间时段。“日”分解模数可以根据数据特点,即工作日和休息日确定为2。“区域”分解模数即为不同类别区域数量k,如果k过小致使不同类别区域冗杂在一起,如果k过大则会使同一类别区域分在不同的模式中;在同一组参数下 SN-TUCKER分解结果歧义性减少,但是仍然具有不唯一性。综上所述本实验对X_O与X_D选择k[3,18]分别进行五组张量分解实验,结果如图6所示,在[3,8]区间中,随着k值的增大,图6（a）和图6（b）中的目标函数呈下降趋势,在[8,18]区间中,目标函数趋于稳定,因此本实验选取k=8,此时目标函数基本达到[3,18]区间的最小值,且具有较强的可解释性,综上所述,本文的张量分解模数确定为8×4×2。经过实验发现,当对A和B进行稀疏约束,不对G和C进行稀疏约束,可获取张量分解理想结果。

显示原图| 下载原图ZIP| 生成PPT

图6   目标函数随着k值的变化

Fig. 6   The change of objective function with k

4.3 张量分解结果分析

本文使用出租车OD数据分别构造上车张量X_O和下车张量X_D,通过SN-TUCKER分解,将张量X_O和X_D分别分解为一个核心张量和3个因子矩阵： A∈$\mathbb{R}^{11 663\times 8}$、B∈$\mathbb{R}^{24\times 4}$和C∈$\mathbb{R}^{16\times 2}$的乘积形式,A、B、C可以被认为是每一特征维度的主成分,即空间维度、时段维度与日维度的主成分;因子矩阵的每一列分别代表特征空间的一种模式,即区域模式、时段模式与日模式;A每一行代表某一空间格网与对应的8种区域模式相关程度,B每一行代表每一时段与对应的4种时段模式相关程度,C每一行代表每一天与对应的2种日模式相关程度,模式值越大相关程度越大。下面分别对每个特征模式的结果进行分析。

4.3.1 日模式

日维度的因子矩阵C∈$\mathbb{R}^{16\times 2}$有2种模式,如图7,其中2012-11-03（周六）、2012-11-04（周日）、2012-11-10（周六）和2012-11-11（周日）是休息日。图7（a）中,红线在休息日的模式值接近0,在工作日的模式值是高值且平稳状态,表明红线在工作日特征信息较明显;蓝线在休息日的模式值达到高峰,而在工作日模式值接近0,表明蓝线在休息日特征信息较明显,因此红线是工作日模式,蓝线是休息日模式。同理,图7（b）中,红线是工作日模式,蓝线是休息日模式。

显示原图| 下载原图ZIP| 生成PPT

图7   日模式变化规律

Fig. 7   Day patterns

4.3.2 时段模式

时段维度的因子矩阵B∈$\mathbb{R}^{24\times 4}$具有4种模式,如图8。图8（a）蓝线与图8（b）蓝线在时段7-8达到高峰期,因此图8蓝线对应早高峰模式。图8（a）紫线与图8（b）紫线在时段9-16期间达到高峰期,图 8（a）紫线在时段12存在低谷现象,图8（b）紫线在时段13存在低谷现象,低谷现象是由于人们处于午休时段,出租车需求略微下降,因此图8紫线对应日间模式。图8（a）红线与图8（b）红线在17-21时达到高峰状态,因此图8红线对应晚高峰时段。图 8（a）绿线与8（b）绿线在时段0附近达到高峰,因此图8绿线为夜间模式。

显示原图| 下载原图ZIP| 生成PPT

图8   时段模式变化规律

Fig. 8   Hour patterns

4.3.3 区域模式

张量分解结果中的因子矩阵A∈$\mathbb{R}^{11663\times 8}$,表示每个格网中对应的8种上下车模式,每个模式值表示对应的模式在该格网中的贡献值,该值与格网的上车量成正比,即格网上车量越大,贡献值越大。由于每个网格的上下车量本身存在很大差异,导致无法直接比较不同格网的模式值,因此本实验对因子矩阵A在行方向进行归一化处理,如式（3）所示。

Ap,q'=Ap,q∑pNAp（3）

式中：分子表示A第p行q列元素,分母表示第p行元素之和。

矩阵A_O归一化后,得到矩阵 AO'∈$\mathbb{R}^{11663\times 8}$,表示 8种模式在对应格网中的贡献比例,比例越大,说明对应的模式在该格网中越突出。根据 AO'中每一行的值,可获取每个格网中区域模式组成,部分格网以1种模式为主,部分格网包括多种模式。将属于每种模式的所有格网在地图上显示出来,得到上车点的8种区域模式,如图9,每一区域模式中的上车点具有相似的时间特征。同理,可获取下车点的8种区域模式,如图10所示。

显示原图| 下载原图ZIP| 生成PPT

图9   上车点区域模式分布

Fig. 9   Spatial patterns of Os

显示原图| 下载原图ZIP| 生成PPT

图10   下车点区域模式分布

Fig. 10   Spatial patterns of Ds

4.3.4 核心张量分析

通过张量分解得到的上下车时间模式,可以判断出租车用户的行为时间特征。结合核心张量,建立时间模式和区域模式的对应关系,可初步推断出租车用户的行为空间功能特征。核心张量Go∈$\mathbb{R}^{8\times 4\times 2}$与G_D∈$\mathbb{R}^{8\times 4\times 2}$,如表2和表3,分别表示8种上车点/下车点模式、4个时段模式和2个日模式之间的联系强度。根据表2与表3可确定上车点和下车点之间的相互作用,建立同一个时段中活跃的上车点与下车点的映射关系（表4）,从可而进一步判定上车点所在区域与下车点所在区域对应的功能。

工作日早高峰时段,O₅与D₅较活跃,该时段人们通常从居住区前往工作区,可以初步判断O₅对应居住区,D₅对应工作区。工作日日间时段,O₈、O₂与D₆较活跃,O₂在工作日晚高峰时段也较活跃,O₂、O₈与D₆空间分布广泛,白天人流较大且活动范围较广,推测O₂、O₈与D₆为餐饮购物、交通、景区、教育、医疗等多种功能混合区域。工作日晚高峰时段,O₇、O₂与O₁上车数量较多,D₂、D₁下车数量较多,其中O₇空间分布较广,推测O₇为餐饮购物、交通、景区、教育、医疗等多种功能混合区域;O₁在休息日晚高峰时段也较活跃,初步推测O₁为景区或娱乐区（酒吧、KTV等夜生活场所）;人们晚高峰时段通常是回家或者到达娱乐场所,初步判断D₁、D₂为居住区或娱乐区。工作日夜间时段,O₃、D₇与D₃较活跃,夜间娱乐场所通常是出租车上下车的热点区域,居住地则会出现大量出租车下车现象,推测O₃对应娱乐区,D₃与D₇为居住区或娱乐区。

休息日早高峰时段,O₆与D₄较活跃,休息日早晨人们活动起点通常是居住地,目的地为景区或者机场车站,推测O₆为居住区,D₄对应景区或者机场车站。休息日日间时段,O₄与D₈较活跃,白天人流较大且活动范围较广,推测O₄与D₈为餐饮购物、交通、景区、教育、医疗等多种功能混合区域。休息日晚高峰时段,O₁与D₁较活跃,上文可知O₁为景区,D₁为居住区/娱乐区混合区域。休息日夜间时段,O₃与D₃较活跃,上文可知O₃为娱乐区, D₃为居住区/娱乐区混合区域。

图9与图10中包含8种上车/下车空间模式分布,表明在不同时段,出租车用户到达不同的场所,进行不同的活动,揭示了出租车用户出行功能区的动态变化。

4.4 出租车用户出行功能区识别

通过张量分解模型得出的时空模式可以初步探究出租车用户出行的目的,有些时段的出行行为目的明确,而有些时段则可能有多种出行目的,难以确切识别目的地的区域功能,因而,本文利用POI类型判断区域功能特征,进一步提高出租车用户出行功能区识别的准确性

4.4.1 区域模式功能特征

本文计算每种区域模式（图9与图10）下的POI频数密度值（Frequency Density,FD）,计算公式见式（4）,其中N_{i, j}表示第i种区域模式的第j类POI数量,S_i表示第i种区域模式总面积。

FDi,j=Ni,jSi（4）

为了消除不同POI类别数量级差异影响频数密度,对FD进行Min-max标准化,计算公式见式（5）,其中FD_min与FD_max分别是同一类POI在8种区域模式中最小和最大频数密度。上车点区域模式与下车点区域模式POI归一化频数密度见表5-6。

FDnor(i,j)=FDi,j-FDminFDmax-FDmin（5）

O₁区域分布比例最高的POI是风景名胜,表明该区域主要为景区,排除了上文娱乐区推测;O₂区域中每类POI的FD_nor都很高,具备大量餐饮、购物、休闲、公司、住宿等服务,具备成熟商业区特征,包括多种功能;O₃区域包含大量体育休闲、餐饮、购物、住宿等服务设施,说明该区域为娱乐区;O₄、O₇与O₈具备丰富的餐饮购物、生活服务、体育休闲、住宿、交通设施服务,具备新兴商业区特征,包括各种功能;O₅与O₆区域具有大量住宅POI,且包含丰富的餐饮、生活、医疗、科教文化等基础设施,具有明显居住区特征。

D₁、D₂与D₇区域具有大量住宅、科教文化、交通设施、餐饮购物等服务设施,同时体育休闲服务丰富,推测D₁、D₂与D₇区域为居住区或娱乐区混合区域,与上文推测一致;D₃空间区域住宅服务POI分布比例最高,同时配套一定数量的餐饮、购物、生活、医疗等基础服务,具备居住区的特征,排除了上文娱乐区推测;D₄空间区域主要包括风景名胜,因此D₄主要为景区;D₅空间区域具有丰富的公司企业、餐饮和交通设施等服务,具备工作区的特征;D₆区域中每类POI的FD_nor都很高,具备大量餐饮、购物、休闲、公司、住宿等服务,具备成熟商业区特征,包括多种功能;D₈具备丰富的餐饮购物、生活服务、体育休闲、住宿、交通设施服务,其FD_nor值略低于D₆,具备新兴商业区特征,包括各种功能。

在核心张量分析出租车用户上车点/下车点对应功能区基础上,结合城市POI数据进一步推测功能特征,二者结果大体一致,同时可进一步揭示出不同时段下的出租车用户出行区域功能特征。

4.4.2 功能区识别

通过以上分析,可判断O₅、O₆、D₃为居住区,O₁、D₄为景区,O₃为娱乐区,D₅为工作区,D₁、D₂和D₇为居住区/娱乐区的混合区域,O₂、O₄、O₇、O₈、D₆和D₈为商业区。为进一步探索北京市六环整体功能,将每一个格网对应的O区域比例向量 AO'和D区域比例向量 AD'合为一个向量 AOD'∈$\mathbb{R}^{11663\times 16}$,根据 AOD',我们可以判断每个格网属于的O区域模式与D区域模式,一个格网可能为一种O区域模式或者一种D区域模式,也可能为多种O区域模式和多种D区域模式的混合模式。将每个格网的区域模式组成对应的功能组合在一起,得到格网的功能组合,结果如图11,其中“未识别区”对应出租车上车和下车数量为0的区域,这类区域通常是待开发区或者出租车无法到达的区域。图11中商业区包括A₁北京CBD、A₂王府井、A₃西单等商业区域。工作区包括B₁中关村科技园、B₂亦庄经济开发区等。景区包括C₁故宫、C₂天坛、C₃颐和园等。居住区往往和其它功能区组合在一起,如工作区、景区、商业区、娱乐区,形成混合区域;从图11中可看出,具有居住功能的区域分布空间较广泛,其中居住区和景区混合区域通常是住宅小区和周围公园的组合（D₁、D₂区域等）;居住区和工作区混合区包括E₁磁器口附近、E₂亦庄科技园附近等;居住区和商业区混合区包括F朝阳门附近等;居住区和娱乐区混合区包括G₁北京工人体育场附近、G₂簋街等;居住区、娱乐区和商业区混合区则包括H₁蓝色港湾、H₂南锣鼓巷、H₃五道口等区域。综上,结合区域活跃时间段与对应POI分布而获取的空间区域语义特征与现实生活基本常识认知基本一致,验证了本文分析结果的合理性。

显示原图| 下载原图ZIP| 生成PPT

图11   功能区

Fig. 11   The functional regions

图11是8种上车点区域模式和下车点区域模式对应城市功能区的叠加结果,即所有时间段内的出租车用户出行功能区的叠加。图9与图10中的空间模式对应功能区以单一功能区为主,而图11中的单一功能的格网区域数量远小于多功能的格网区域数量,表明一个区域往往不是静态的功能,而是随着社会人群的需求变化而变化。

5 结论

本文根据北京出租车OD轨迹数据,基于城市人群在不同时间段进行不同活动,使用张量分解模型探究北京城市人群在时间上与空间上的出行规律。实验结果表明,利用张量分解得到的出租车OD数据时空模式,可以从不同维度的特征挖掘出出租车用户出行时空活动模式,从而识别不同时段中的出租车用户出行功能区。本文利用城市人群出行模式获取不同时段的城市功能区分布,并将不同时段的功能区进行叠加,获取所有时段中的综合功能区分布,其混合功能区数量大于每个时段的混合功能区数量,表明城市区域功能不是静止不变的,而是不同时段功能的组合。但是根据出租车OD数据识别的用户出行功能区有限,需结合POI或其它富含语义信息的数据,提高上下车点所在区域功能信息判断的准确性。

城市人群的出行不局限于出租车,因此出租车不能表示所有城市人群的出行信息,同时出租车上下车点局限于路网区域,即不能到达所有研究区域,因此出租车OD数据隐含的语义信息不能全面反映出城市所有功能区的时空分布。为了进一步识别出人群出行功能区,本文计划在出租车轨迹数据基础上,融合其它来源的数据,如社交媒体定位或文本数据,从而更加全面而准确识别城市人群出行功能区分布,挖掘出城市人群出行驱动机制。

The authors have declared that no competing interests exist.

参考文献

文献选项

• 原文顺序
• 文献年度倒序
• 文中引用次数倒序
• 被引期刊影响因子