模式指对象活动的一般方式, 数据中反映的模式通常用于描述用户的行为规律.Zhang等^[7]提取高校学生的行为模式特征, 分析其行为与课程绩效的关联性.在金融领域中, Huang等^[8]针对区块链中模式, 提出相似性度量, 用于高精度聚类.在Web应用中, 模式用于发现用户的特征、偏好和习惯, Chovanak等^[9]将典型重复的用户行为特征作为行为模式, 实现个性化预测和推荐等任务.同样在Web浏览行为分析上, Hussain等^[10]利用用户上网行为模式和网络轨迹, 分析用户兴趣和认知预测.

时空轨迹数据记录用户的位置序列, 近年来学者考虑将出行路径相似的用户进行聚类, 再发现其共同的行为规律.部分研究利用主成分分析提取轨迹数据中稳定的移动模式, 用于用户位置预测或是静态POI对GPS轨迹数据的标注, 再进行语义判断^{[11, 12]}.Cao等^[13]提出基于表示学习的用户轨迹语义表示系统, 通过向量的形式表示轨迹数据特征.Maiti等^[14]针对空间数据中提取的行为模式, 提出相似性度量及对应模式的挖掘方法.

对于多源异构的轨迹数据, 经过融合集成后往往存在高维度和高稀疏性的问题, 无法直接进行模式抽取分析, 通常需要进行降维.Xu等^[15]总结大规模异构数据处理时经典的基于机器学习的数据降维方法及其应用.Yan等^[16]开发面向轨迹文本数据的特征分析算法, 用于提取数据中的特征, 达到有监督或无监督的数据降维目的.

城市功能区划分是通过POI分布或用户行为发现具有某种主题功能分布的空间区域.随着数据获取途径的增加及GPS的普及, 利用用户轨迹数据进行功能区划分成为主要方法.功能区划分方法主要分为基于统计学划分和基于机器学习划分两种.

基于传统统计学方法, Hu等^[4]提出利用城市POI数据的分布和其功能类型比值等统计特征划分城市功能区, 并分析不同功能区的空间分布特征.Qi等^[3]根据出租车乘客上下车的时间变化规律进行区域功能判断, 使用聚类的方法划分功能区域.

机器学习技术也被应用于城市分析和功能区划分.Tao等^[17]利用GPS数据, 构建四维张量, 发现用户行为在不同时间的规律, 并进行空间聚类以划分城市区域.Papadakis等^{[18, 19]}量化空间信息的语义, 对地理空间赋予功能属性, 整合至地理信息系统中, 进而结合POI属性和数据驱动, 提出基于隐含狄利克雷分布(Latent Dirichlet Allocation, LDA)的功能区划分模型, 并评价划分结果.在此基础上, Wang等^[6]为了理解多中心城市区域的空间结构, 提出基于狄利克雷多项式回归的模型, 利用交通模式和POI数据获取城市单元功能, 发现不同区域之间的相互关系及城市的功能分布.Yuan等^[1]结合POI数据和出租车轨迹数据, 基于主题的LDA模型提取用户出行规律, 并分析各功能在不同地理位置的强度.

总之, 目前研究大多从空间维度上利用不同方法进行城市功能的划分并得到功能区.

为了从时空数据集中获取有价值的信息, 时空聚类能划分具有相似时间、空间的数据集合, 发现时空要素的发展规律和趋势.

Parente等^[20]旨在将具有相似属性的对象分组到各自的类别中, 提出时空排列扫描统计方法, 完成时空聚类任务.Paci等^[21]提出以状态空间为中心的贝叶斯半参数混合模型, 完成时间序列的聚类, 引入时空变化的混合权值, 设计基于时空模型的方法, 实现时空数据中的动态聚类, 预测时空依赖性.Meliker等^[22]利用依赖于时间最近邻的方法, 对用户居住史进行时空聚类分析, 发现不同类别之间的潜在关系.

许多时空现象具有顺序性和层次性, 研究对象具备固定的空间位置和规则的时间步长.传统的时空周期模式挖掘并未考虑其内在的层次性, Zhang等^[23]提出基于Traclus的轨迹聚类算法, 补充时空序列, 在轨迹数据聚类时同时考虑方向、速度等语义时空信息.Pereira等^[24]提出时空排列扫描统计量, 识别聚类结果的数量、位置和时空维度, 提高聚类结果的精确度和鲁棒性.

本文以北京市主城区作为分析区域, 网格化处理空间区域, 得到多个相同规格的矩形空间单元.对用户轨迹数据融合并标注北京POI访问次数分布情况, 标注的方法为基于概率传播的匹配模型^{[25, 26]}.统计空间单元中各功能类别POI访问次数后获得单元功能主题特征, 其中特征向量的一个维度表示一个主题功能.

由于POI类别的多样性和轨迹数据的稀疏性, 单元功能主题特征分布无法准确表示单元空间内用户的活动规律.为了解决此问题, 本文基于自编码器的表示学习获取数据的潜在低维特征, 用于准确发现单元内部用户活动规律.

自动编码器神经网络是一种应用反向传播的无监督学习算法^{[27, 28]}, 对于一组未标记的训练样本集V_n={x₁, x₂, …, x_n}, 自编码器不断学习函数h, 使得h_{w, b}(x)≈ x, 调整权重项W和偏置项b, 使目标输出值等同输入, 即y_i= x_i.本文使用增加稀疏惩罚项的自编码器, 目的是找到隐藏层中神经元集合s数据内部的规律.方法需要对s中部分神经元进行抑制, a_j(x_i)表示隐藏层中当输入x_i时神经元j的激活度, 可计算该神经元在不同输入下的平均活跃度:

φ︿j= 1n∑i=1na_j(x_i).

为了保证 φ︿j=φ , 将两者之间的KL散度(Kullback-Leibler Divergence)作为正则项, 加至原损失函数中进行稀疏限制, 通过系数β 调整限制幅度, 修改后的损失函数为:

J_sparse(W, b)=J(W, b)+β ∑j=1sKL( φ︿j‖ φ ).

稀疏自编码器可有效表示高维稀疏的数据, 同时能提升线性可分性.

基于自编码器的降维模型输入为空间G内标注后的高维行为数据.为了解决轨迹数据稀疏的问题, 本文选择在模型输入时考虑空间邻域信息, 使单元空间向外扩展, 并将其映射至同一矢量空间进行空间维数压缩并表示.由于城市功能的高度聚集性, 对于扩大范围后的数据, 降维模型能有效降低数据稀疏和噪声数据的影响.整个学习过程如图2所示.

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 空间单元特征数据嵌入表示学习过程Fig.2 Representation learning process of embedding feature data of spatial units

对于空间单元g_i, 将其空间邻域网格区域g_n(n=1, 2, …, e)数据同g_i的主题特征合并为ρ +e× ρ 维向量进行输入.

观察发现, 同空间区域内, 部分相邻时间跨度的主题特征具有相似分布, 需要将向量的细粒度划分合并为粗粒度划分, 用于发现行为规律.本文提出基于嵌入表示的功能主题模式抽取算法(Thematic Patterns Extraction Based on Embedding, TP2E), 将单位时间跨度聚类后对整体区域进行空间划分, 完成时空聚类后, 得到功能主题模式.算法1描述功能主题模式抽取算法的过程, 分为时间维度聚类和空间维度扩展两部分内容.

算法1 TP2E

输入 时段T内的空间G={g₁, g₂, …, g_n}

输出 功能区主题模式集P

1.初始化F={ }, C={ }, P={ };

2.将g_i同邻域空间合并后进行表示学习;

3.for (每个空间单元g_i)

4. 获取g_i在T内功能主题特征向量表示f;

5. F_i=Temporal_Cluster(g_i, f ) ;

6. 将主题特征矩阵F_i加入集合F中;

7.end for

8.基于F中时间跨度划分对空间单元g_i聚类, 得到

簇C={c₁, c₂, …, c_m};

9.for (每个簇c_i)

10. P(R, T')=Spatial_Extension(F, c_i) ;

11. P=P∪ P(R, T');

12.end for

13.return P

算法1的执行过程如下.第1~2行将输入的时间段T内最小空间单元g_i的主题特征向量进行嵌入表示.第3~7行调用算法2实现时间维度的聚类, 获得重划分时段T内g_i特征矩阵F_i并加入集合F中.第8~13行对所有空间单元, 基于时段划分, 利用K-means聚类后得到簇集合C, 调用算法3在每个簇内对g_i实现空间扩展, 返回区域划分R和对应时段划分T'的集合, 计算集合内对应的主题分布向量, 并入结果集P中.

K-means效果较优, 收敛速度较快, 对于时段划分数据较契合, 能通过误差平方和(Sum of Squares due to Error, SSE)肘方法确定K值.

算法2为一个时间维度聚类算法, 对每个空间单元上的主题特征向量进行聚类.

算法2 Temporal_Cluster(g_i, f )

1.while(∀ 相邻时间跨度向量JS(f_t, ft+1)≤ θ )

2. 合并向量并更新原本跨度为t'=[t, t+1];

3.得到g_i在时间段T上的跨度重划分;

4.计算g_i的主题特征向量矩阵F_i;

5.return F_i

算法2执行过程如下.第1~2行判断该单元内在时间跨度t内的主题向量f_t和其相邻时间点的向量 ft+1是否相似, 若相似合并两者, 更新时间跨度为t', 迭代合并, 直至空间g_i上不存在相似主题特征向量.将合并后的主题特征矩阵F_i返回至算法1中.

在空间单元g_i上的部分主题特征向量会发生时间演变.若相邻时间跨度的两个特征向量相似, 呈稳定状态, 需要对向量平均化处理并更新时间跨度.否则表示主题功能结构存在较大差异, 不需要合并.算法2最终会输出各空间单元在新划分时间段T内的低维特征矩阵, 用于后续进行同时考虑时空维度上的区域空间扩展.

本文选择基于KL散度演变的JS散度度量时间演变后的两个功能主题特征之间的相似度, 见算法2第1行.JS散度定义如下:

其中p_i(i=1, 2)表示d维的概率分布向量.

JS散度结果为[0, 1]的常数, 对于2个向量完全没有重叠时的结果仍有意义, 并解决KL散度非对称的问题, 有效应对稀疏数据中概率为0的情况.

算法3为功能主题模式抽取任务中的空间扩展过程.本文提出空间扩展算法, 进行空间区域划分.该算法在具有相同时段划分的簇内寻找最相似的2个相邻空间单元, 记为初始空间N, 如图3(a)所示.阴影区域作为候选区域, 基于相似度判断是否合并, 不断进行空间边界扩张, 直至簇内剩余空间单元无法合并, 且个数小于给定可忽略值γ (如图3(b)所示), 完成扩展, 得到时空划分集合.可忽略值γ 为簇内数据量过小的单元个数, 取值区间为[0, c_num/2) , 其中c_num表示簇内单元总数.算法3具体步骤如下.

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 空间扩展算法过程Fig.3 Process of space expansion algorithm

算法3 Spatial_Extension(F, c)

1.初始化R={ }, T'={ };

2.while (c≠ Ø ) do

3. 计算簇c内相邻空间单元相似度;

4. for (c内空间单元数可忽略值大于等于γ ) do

5. 合并簇c内最相似空间单元, 记为N;

6. 空间N的横、纵相邻区域加入候选空间集合N_c中;

7. if (Optimized_Merge(N_c, N)≤ θ ) then

8. N同最相似候选空间合并后更新;

9. 从c内移除候选空间单元信息;

10. else

11. 将空间N存入集合R中;

12. 将所在簇的时段划分存入集合T'中;

13. end if

14. end for

15.return(R, T')

算法3执行过程如下.第4~5行寻找簇内某两个最相似空间, 合并主题特征矩阵后记该空间区域为初始空间N, 利用最小距离法确定相似度阈值θ .最小距离法为分类器内确定分类阈值的基本方法, 求解未知类别向量q到事先已知的各类别中心向量的距离以确定阈值θ .第6行将初始空间N的矩形邻域作为候选空间并存入集合N_C中.第7~13行将候选空间集合同空间N进行相似度计算, 选择最相似的两个空间合并, 更新空间N作为合并结果, 移除簇内的空间信息, 完成区域扩展.最后将对应空间划分和时段划分存入集合R和T'中.

对于同簇内不同空间单元来说, 任意2个单元时间段划分仍有不同, 为了更准确地对这种空间单元进行相似性度量, 本文提出优化合并度量的计算方法, 见算法4.

算法4 Optimized_Merge(g_x, g_y)

1.获取两个子空间的主题特征矩阵F_x、F_y;

2.在各时间跨度内随机取时间点, 加入临时序列L中;

3.按照序列L中的时间点从F_x、F_y抽取向量, 得到临时对比矩阵F'_x、F'_y;

4.基于F'_x、F'_y计算相似度ω ;

5.return ω

算法4执行过程如下.第1~2行对于两个相邻空间单元g_x、 g_y, 根据其不同的时段划分, 在各时间跨度内随机取点, 记为临时序列L.第3~4行对于这两个相邻空间分别按照L中的时间点取出对应向量, 得到两个规格相同的概率分布矩阵后, 计算余弦相似度.算法4能将具有不同时间跨度的空间进行向量抽样, 基于抽样结果进行相似度对比.

相比现有模式发现算法, TP2E考虑模式时间演变问题.获取特征矩阵F后, 将n个空间单元聚类为k个簇, 实现时空聚类的代价为

|F|× nk× nk,

可表示为O(|F|n²), | F|为特征矩阵F的大小.

实验数据包括地理空间数据集和用户访问轨迹数据集.地理空间数据集为北京市187万条POI数据, 包含地理坐标、地址、所在区域等属性.用户访问轨迹数据来自百度地图查询数据和摩拜共享单车轨迹数据这2个真实数据集.

百度地图查询与共享单车数据集上包含用户出发时间、起始位置和终点位置等属性, 数据集样例如表1所示, 起始位置由经纬度表示.百度地图查询数据为2017年5月北京市数据, 本文实验所选区域涉及用户数约72万个, 数据量为1 243万条.摩拜共享单车数据为2017年5月北京市单车骑行轨迹数据, 实验选择区域涉及用户数约61万个, 数据量为321万条.

表1

Table 1

表1(Table 1)

表1 实验数据集样例 Table 1 Examples of experimental datasets 数据集 出发时间 起始位置 终点位置 百度地图 查询数据 2017年5月20日 星期五 11:08:12 116.351781 39.944615 116.454307 39.773794 摩拜单车 轨迹数据 2017年5月14日 星期六 15:22:30 116.092269 39.941483 116.27586 39.912643

表1 实验数据集样例 Table 1 Examples of experimental datasets

本文选择如下功能区划分模型进行实验对比.

1)BaseDC.基于密度的功能区划分算法, 确定区域中的子空间后, 计算空间内功能频率密度和功能类别内外排序, 通过POI密度数据发现城市功能分布特征^{[2, 4]}.

2)BaseLDA.基于LDA的功能区划分算法, 在利用LDA提取城市功能结构时, 将空间区域视作文档, 将区域内功能视为主题, 而对应POI类别视为元数据, 方法类似于文本上下文中的主题抽取.形成三层贝叶斯模型架构后, 基于LDA模型能使用一个概率分布表示该区域的功能结构, 发现主题关联, 实现功能区划分^{[1, 6, 19]}.

3)TPE.本文提出的无降维的功能主题模式抽取算法, 对空间单元数据实现时空聚类后获取主题功能模式, 基于主题模式获取功能区划分.

4)TP2E.本文的功能主题模式抽取算法, 在完成对空间单元数据降维后进行时空聚类, 对各空间单元的时间跨度进行重划分后实现空间扩展, 发现具备时间敏感的功能主题模式, 获取动态城市功能区域的划分.

本文选取北京市朝阳区作为实验区域, 空间面积为470.8 km², 需要对空间网格化处理.现有相关研究将区域划分为0.5 km× 0.5 km至1 km× 1 km的规则网格^{[2, 4, 17]}, 本文选择将实验区域划分为1 km× 1 km的规则网格, 单位时间跨度t为1 h, 选取部分区域展示模型有效性.

图4为数据嵌入模型的损失曲线变化, 自编码器降维模型可有效降维并解决数据稀疏问题.

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 自编码器降维模型的训练效果Fig.4 Training effect of dimension reduction model for AE

图4(a)中无处理数据在自编码器降维模型中的损失呈曲线重叠状态, 对于均匀分布的稀疏数据, 验证损失往往高于训练损失, 而且使用数据量稀疏的空间单元发现功能聚集区域的效果较差, 故需要合并单元的空间邻域.将空间单元与向外一周所有邻域合并, 作为输入进行降维, 损失(如图4(b)所示)趋于收敛, 模型的学习效果较优.选择向外两周的邻域空间单元同目标空间单元合并, 损失(如图4(c)所示)存在波动, 故选择第1种数据嵌入表示方法.

基于主题模式的功能区划分可归属为时空聚类任务, 因此可利用聚类性能评价指标进行评价.本文选择两种典型的聚类评价指标— — 轮廓系数(Silhouette Coefficient, SC)和Calinski-Harabasz系数(CH), 并计算SSE, 利用SSE肘方法发现最优时间跨度聚类个数k, 完成K-means聚类任务.

数据嵌入表示在解决数据稀疏的同时提升聚类性能, 故TP2E聚类效果具有明显优势.不同聚类数量k下的度量值变化如图5所示.

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 簇个数k变化时各模型的聚类性能对比Fig.5 Clustering performance comparison of models with different numbers of cluster k

各模型的最优聚类数量k及对应指标值如表2所示.由表可见, 基线方法未能解决数据稀疏问题, 导致无法精准选取特征, 不同簇间相似度较高.TP2E增强空间单位对邻域的相关性, 解决数据稀疏问题后, 聚类效果具有明显提升.

表2

Table 2

表2(Table 2)

表2 各模型的最优k值及指标值 Table 2 Optimal k value and index value of different models 模型 CH系数 k SC系数 k BaseDC 17.92158 5 0.204015 7 BaseLDA 16.36444 3 0.200354 3 TPE 17.865752 3 0.230989 5 TP2E 29.545715 6 0.479502 5

表2 各模型的最优k值及指标值 Table 2 Optimal k value and index value of different models

不同周期的相同时段数据合并后的SSE值如图6所示.针对SSE曲线, 通过肘方法取SSE曲线中下降率变缓时的k值, 判断最优时间跨度聚类个数, 并通过CH系数和SC系数验证后, 按照该值对局部区域的时间进行聚类, 获得不同时间跨度.

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 k不同时的SSE值Fig.6 SSE with different k

本文选择功能频率密度和比率指数对划分城市空间区进行可视化, 用于评价划分结果.频率密度

FD_i= aiAi,

其中, i=1, 2, …, ρ , 表示功能类型, a_i表示时空单位内i类型的数据量, A_i表示时间单位内i类型的总数.

比率指数

RA_i= FDi∑i=1ρFDi× 100%, i=1, 2, …, ρ ,

为类别为i的数据在空间单位内的频率密度之比.不同空间内各功能比率指数RA≤ 40%时, 认为是混合功能类区域.本文利用3.2节进行功能主题模式的抽取, 计算各模式的比率指数, 确定各区域的功能, 并以图形化的方式展示.

各模型的功能区划分结果如图7所示, 图中颜色区分不同主题功能.基于主题模式划分的功能区在不同时间跨度具有明显差异.基于密度的功能区域划分结果如图7(a)所示, 在功能空间分布上较分散, 无法获得合适的城市功能区域.基于LDA的功能区域划分结果如图7(b)所示, 能划分合适功能区域, 但对于整个时间周期的区域功能划分粒度过粗.

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 各模型的功能区划分结果Fig.7 Results of functional zoning by different models

TP2E获得的9时~10时和12时~15时这2个时间跨度内的功能区划分结果如图7(c)、(d)所示, 在满足一定划分粒度的前提下, 仍能划分较真实的功能区域.对于同一区域, 图7(c)中可看出餐饮类和购物类功能区占比最高, 而在图7(d)中可发现在其它时间跨度变为购物类功能区占主导地位, 即发生城市功能结构时间演变.这种演变是无法通过现有功能区划分方法体现的, 而利用功能主题模式发现的功能区划分方法能更好地诠释这种变化规律.