一种聚类隐马尔可夫模型的时空轨迹预测算法

孙 红，陈 锁

1(上海理工大学，上海 200093) 2(上海现代光学系统重点实验室，上海 200093)

1 引 言

近年来，由于大数据时代的来临，以及物联网技术的广泛应用，各种移动终端的大量普及，大量GPS数据随之产生，各种轨迹数据也越来越容易获得，对GPS轨迹数据的研究也越来越多.大部分终端设备都已装有GPS模块，这些终端设备因而可以记录GPS位置序列，这些序列包含经纬度，时间戳等信息.时空轨迹(Trajectory)是移动对象的位置和时间的记录序列.作为一种重要的时空对象数据类型和信息源，时空轨迹的应用范围涵盖了人类行为、交通物流、应急疏散管理、动物习性和市场营销等诸多方面.时空数据的急剧增加，加之时空数据处理更为复杂，使数据处理任务日趋繁重的形势更加严峻.因此，寻找有效的时空数据处理方法具有十分重要的意义.

时空大数据与非空间数据相比，具有空间性、时间性、多维性、海量性、复杂性等特点.关于时空大数据，目前存在一些技术前沿领域，如云计算方法和挖掘技术.而预测移动对象的运动趋势有着广泛的应用.比如能够有效地在轨迹数据中发现隐藏的有价值的信息，如频繁路径、个性化兴趣点以及移动意图等.这些结果将对未来的商业应用以及管理活动提供有力的信息保障.

本文通过分析时空轨迹数据，提出一种根据移动对象经过的轨迹，预测移动对象下一个轨迹点，以及未来的终点位置的模型.

2 研究背景

对于移动对象轨迹位置预测，近年来主要产生了两种类型.一种是对历史轨迹频繁模式进行挖掘进而进行预测.另一种是基于马尔可夫(Markov)模型，对历史轨迹进行建模并预测.

轨迹频繁模式挖掘算法，目前主要分为两种：Apriori算法和PrefixSpan算法.一些研究将它们应用到轨迹位置预测的建模中.一种基于广度优先搜索的Apriori算法结合自定义的运动函数和移动轨迹模式，堪称一种创新型的方法，可以参考Jeung等人[10]的相关文章.还有一种改进版的Apriori算法，可参考Morzy[11]的相关文章.另一种主流算法是PrefixSpan算法，文献[12]提出一种改进型的PrefixSpan用来发现运动轨迹的历史频繁项并生成预测模型.此外，Monreale等人[2]提出的WhereNext方法也不失为一种好的方法，该算法从历史轨迹中提取具有几种不同运动行为的运动模式，进而生成预测模型.

马尔可夫模型对轨迹序列建模中，有一些比较新的成果，比如Figueiredo等人[7]，他们使用不稳定，瞬态且时间多样化的轨迹序列，构造出一种非参数模型用来对轨迹位置等做出预测.文献[9]使用了一种HITS-based模型来挖掘用户轨迹兴趣点，侧重于分析轨迹全局的运动模式.

总的来说，现有工作大多是依据大量移动对象轨迹的重合部分，挖掘出轨迹相似行为，再用来做预测.这种方式的缺点是：着重考虑表面发生的行为，对于事物内在的隐含信息的考虑有所欠缺.

由于移动物体产生的轨迹数据量比较大，而且这样的数据还在不断地增长，因而对于以上两种方法，模型的训练成本非常大.本文基于此原因，提出一种基于聚类分区域的隐马尔可夫模型，经过大量GPS数据的训练，得出最佳模型参数，最后通过隐状态之间的转移矩阵预测移动对象轨迹最终位置.

3 聚类分区域的隐马尔可夫模型

3.1 隐马尔可夫模型(HMM)

马尔可夫模型中，状态是可见的，因此其又被称为可视马尔可夫模型(visible Markov model，VMM)，状态可见，在一定程度上限制了模型的适应性.而在隐马尔可夫模型(HMM)中，每个状态是不可见的，可见的是各个状态产生的事件的概率分布，状态之间的转换也是不可观察的，可观察的事件是由状态产生的，每个状态产生一个特殊分布的事件集合.图1描述了隐马尔可夫模型的基本原理.

图1 隐马尔可夫模型Fig.1 Hidden Markov model structure

一个HMM由如下几个部分组成：

1)模型中隐状态的数目N；

2)隐状态可能输出的不同可观察符号的数目M

3)状态转移矩阵A={aij}，其中，

aij=P(qt=sj|qt-1=si)，1≤i，j≤N

aij≥0

4)可观测事件概率分布矩阵B={bj(k)}，其中，

bj(k)=P(Ot=vk|qt=sj)，1≤j≤N；1≤k≤M

bj(k)≥0

5)初始状态概率分布π={πi}，其中，

πi=P(q1=si)，1≤i≤N

πi≥0

一个HMM一般来说可以用以一个三元组来描述，μ=(A，B，π)，A为状态转移概率，B为符号发射概率，π为初始状态概率分布.当一些现象产生是由于某些事物本质变化做决定时，用HMM进行建模一般是非常有用的.因为HMM能够“深刻理解”事物的变化规律.HMM三个基本问题及其解决方法如下：

1)给定一个模型μ=(A，B，π)，如何快速计算出观察序列O=O1O2…OT的概率，即P(O|μ).

解决方法：前向算法、后向算法.

2)给定一个模型μ=(A，B，π)，和一个已经产生的观察序列O=O1O2…OT，如何找出最可能产生观察序列O的状态序列Q=q1q2…qT.

解决方法：维特比算法.

3)给出大量观察序列集，如何构造出一个合适的隐马尔可夫模型μ=(A，B，π).解决方法：Baum_Welch算法

HMM在自然语言处理研究中有着非常广泛的应用，尤其是在语音识别领域，一段语音可以作为观测序列，语音对应的文字作为隐状态，语音识别的首要任务就是求出最有可能发出该语音序列的隐状态，将语音恢复成文字.通过一部分可观测序列，可以预测对应的隐状态序列，基于此，本文将HMM模型移植到时空轨迹序列预测的问题中.将大样本时空轨迹序列，作为HMM的可观测序列，将局部的时空坐标点看成是附近某个隐状态所发射出的观测值，该隐状态可以视为该局部区域的地标点.如图2所示.

图2 隐状态及观测值Fig.2 Hidden state and visible state

图中黑色三角形即为隐状态，表示城市内的多个地标，这些地标可以看成是人群最爱去的地点的集合.以A点为例，周围有1，2，3，4，5，6几个GPS坐标点，它们分别是四条GPS轨迹中的点，因为都处于A点附近，因此可认为它们是A状态发射出的观测点，即它们服从A的发射概率分布.本文的任务之一就是要建立这样的HMM模型，模型训练完成之后，输入要预测终点的一段时空轨迹序列，即为观测序列，然后通过相应的算法求解出最佳隐状态序列，通过隐状态序列最后一个状态以及转移矩阵，计算出概率最大的下一个隐状态，再根据该隐状态的发射概率分布计算出最有可能的发射值即为预测的终点.整个流程如图3所示.

3.2 聚类分区域的隐马尔可夫模型——CSHMM

HMM模型常用来进行语音识别，一段语音的音节数量是很少的，一般是几十个，这样的数量级，采用维特比算法来求解隐状态，即推断语音对应的文字或单词，算法的执行时间是非常快速的.其原因是维特比算法的时间复杂度是O(N2T)，其中N为发射出来的可观测点的数目，T为状态序列的长度.由于车辆一般都是往一个大致的方向行驶(一般不会存在绕着整个城市或局部区域反复行驶)，因而一条完整的行驶路径所穿过的地标区域个数是比较少的，因而用维特比算法也是很快可以计算出最佳隐状态序列.然而整个城市的所有地标区域是非常多的，比如商场，运动场，小公园，小吃街等，实际上任何地方都可能成为车辆停止的终点.因为时空序列终点的预测精度只要精确到一定的范围内即可，因此所有隐状态的数量就非常多了，可能会有几万甚至几十万个.如果针对整个城市建立一个含有如此庞大数量隐状态的隐马尔可夫模型(尤其是北上广深这样的特大城市)，那么在用维特比算法求解隐状态序列时，递归算法将会非常耗时，甚至很可能导致程序直接奔溃.本文针对该情形，设计了一种聚类分区域的隐马尔可夫模型(Clustering Subregion Hidden Markov Model，即CSHMM).

图3 算法流程图Fig.3 Algorithm flow chart

3.2.1 产生隐状态集

建模的第一步，是要确定所有的隐状态集合，在本文中即为确定整个区域中的地标点集合.

一般来说GPS经纬度坐标数值精确到小数点后4位即可明显的区分两个坐标点为两个不同的地点(相差几十米至150米的两个地点本文当作同一地点)，统计训练数据集中的所有不同地点出现的频率fi.设置一个阈值r，将fi>r的所有地点标记为地标.每个地标使用其频率作为权重.

其次，将得到的所有地标点集，根据其权重进行聚类，聚类的原则是：权重越高的地段所包含的区域面积越小，以使多个子区域包含的数量级相当.在此我们采用坐标点的数量来代替区域面积，可以对此建立一个线性模型(也可以是非线性模型)来确定每个区域内的大致坐标点数目.本文采用一种改进型的kmean——kmeans_cs如下函数模型来判定每个地标周围的坐标点是否属于该地标的周边区域：

设o为一个聚类中心点，p为需要归类的点，n为o类样本集中包含p点的所有轨迹数量，o和p之间的距离D(o，p)采用如下定义：

D(o，q)=e‖o，p‖·‖o，p‖

(1)

其中||o，p||为o，p两点的欧氏距离.由公式(1)可以看出两个坐标点的距离越近，则会有一个大于1的权重乘进去使其增大.如果距离越大，那么乘上的权重会以指数增大，使得距离增加的更快.这样会使得聚类的结果中面积越小的子区域包含的点越多，面积大的区域包含的点越少.改进后的聚类算法可参考下文中2.2.2节中的聚类算法.

图4大致描绘出映射建立结束后的图形.

图4 分区域后的隐状态Fig.4 Hidden state after subregion

图5中的黑色三角形即为隐状态，周围的点即为隐状态的发射符号集.

直观表示如图5所示.

图5 隐状态和观测值的映射Fig.5 Mapping of hidden states and observation values

如图4所示，首先通过聚类算法将城市中的所有地标点进行聚类，以此将全局区域分成若干个小区域，然后针对每个小区域，进行HMM建模.模型建立完成以后，根据输入的待预测的部分GPS轨迹，首先进行区域匹配，选出最佳HMM模型，然后使用维特比算法针对该模型进行预测.

3.2.2 聚类划分子区域

由于地标点数目太多，因而我们需要划分子区域，使得每个子区域内的地标点根据其权重均匀分布.划分的原则是：子区域之间GPS轨迹跨越尽可能少.本文针对GPS轨迹的变化规律设计了一种改进的k-means聚类算法，传统的k-means算法采用欧氏距离来判断点和点之间的类似程度，此方式会打破GPS轨迹坐标的连贯性，即会大量出现同一条轨迹存在于多个子区域中，这是不符合要求的，因为本文需要在每个子区域单独求解HMM模型，每个模型模拟了一个特定子区域GPS轨迹模式.基于该目的，本文提出一种能够按照GPS轨迹来聚类的算法，该算法改进了点和点之间的距离公式：

(2)

设o为一个聚类中心点，p为需要归类的点，n为o类样本集中包含p点的所有轨迹数量，o和p之间的距离D(o，p)采用如下定义：

D(o，p)=tanh(n)×‖o，p‖

(3)

其中‖o，p‖为o，p两点的欧氏距离.采用tanh函数来作为权重函数.不采用sigmoid是因为sigmoid自变量取值为0时函数值为0.5，而这里想要的是在区域p内含有q点的轨迹数目为0时，权重就变为0.

改进后的聚类算法如下：

1)任选K个初始聚类中心Z1(1)，Z2(1)，Z3(1)，…Zk(1)；

根据样本间的相互距离，将样本分配到K个聚类区域，即若Min{d(Zi(k)，X)，i=1，2，…K} = d(Zj(k)，X) = Dj(k)，则X∈Sj(k)，其中，k为迭代序号.

2)更新各个聚类中心Zj(k+1)，其中j=1，2，3...K；

3)若Zj(k+1)≠Zj(k)，则跳转到步骤(2)，将各个样本重新分类，如此迭代；如果Zj(k+1)=Zj(k)，j=1，2，…K，算法收敛，计算完毕.

3.3 分区域进行HMM的模型训练

对于每个子区域，其所有的地标数据已经确定，每个地标周围的称为发射符号集的GPS坐标点也已经确定.已知的信息量比较充足，比如可以根据每个GPS坐标的出现频率，计算出发射概率矩阵B0.也可以根据大量轨迹来计算地标与地标之间的转移概率，从而计算出转移矩阵A0.

Baum-Welch算法是一种EM算法，专门用来训练HMM模型.EM过程保证算法一定能够收敛到一个局部最优点，但是它不能保证找到全局的最优点.所以Baum-Welch算法不能保证找到全局的最优模型.但是由于我们已经知道了大致的解：状态转移矩阵A0和发射概率矩阵B0，如果将它们作为迭代的初始值，那么会有很大的概率收敛到全局最优点，并且训练的时间也大大减小.本文采用的Baum-Welch算法.

对于给定的隐马尔可夫模型(HMM)μ和已知的观察序列O=O1O2…OT，假设在t时刻该模型位于状态si，在时间t+1时刻模型位于状态sj，此概率ξt(i，j)可由如下公式计算获得：

(4)

给定HMM的参数μ和观察序列O=O1O2…OT，在时间t位于状态si的概率γt(i)有下面的公式计算获得：

(5)

而μ的参数可以有下面的公式重新估计：

(6)

(7)

(8)

步骤1. 随机给参数πi，aij，bj(k)赋值，并使其满足如下约束：

由此得到新的模型μ0，令i=0，并用期望最大化算法(EM)重新估计模型的各个参数；

步骤2. EM计算：

E步：由模型μi根据公式(4)和公式(5)计算出ξt(i，j)和γt(i)；

M步：用E步得到的期望值，根据公式(6)，公式(7)，公式(8)重新估计参数πi，aij，bj(k)，得到模型μi+1；

步骤3. 循环计算：

令i=i+1.重复执行EM计算，直到πi，aij，bj(k)收敛.

3.4 使用模型进行预测

对于给定的部分GPS轨迹序列，根据最后一个坐标的位置，确定相应的子区域，选取相应子区域的HMM模型.将不在该子区域内的坐标点切除.使用靠近尾部的坐标点作为可观测序列.问题就变成了HMM模型的第二个问题.因此本文采用维特比算法来求隐状态序列.算法如下：

步骤1. 初始化：

δ1(i)=πibi(O1)，1≤i≤Nφ1(i)=0

步骤2. 归纳计算：

记忆回退路径：

步骤3. 终结：

步骤4. 路径回溯：

求出隐状态序列后，根据转移矩阵A=[aij]和最后一个隐状态qlast，求出下一个概率最大的可观测坐标点Pnext.求法如下：

(9)

其中last为当前时刻，next为下一时刻.

4 实验与分析

本文使用的数据集为T-Drive Trajectory data sample，来自微软亚洲研究院(MSR)，一共有10357个北京出租车在一周内的GPS轨迹，每条轨迹提取出包含经纬度和时间戳信息的子信息.实验选区75%的数据作为训练数据集，25%的数据作为测试数据集.代码使用python语言，硬件配置为 Intel(R) Core(TM) i5-3320 2.60GHZ，12.0GB内存，320GB硬盘的计算机.

实验结果采用准确率，召回率为评估标准.

4.1 本文模型使用kmeans和改进型kmeans结果对比

结果如图6所示(准确率和召回率为改进kmeans后的结果，准确率2和召回率2为不使用改进的kmeanss).

图6 kmeans和改进型kmeans对模型效果的影响Fig.6 Effect of kmeans and improved kmeans on the model effect

结果分析：改进kmeans后的算法导致模型准确率和召回率提高了15%左右，说明改进后的kmeans认为同一个轨迹中的点更应该属于同一个聚类，因为这些点之间是有关联的，有关联的点之间的跳转，直接导致转移概率的增大，因而会有更优秀的准确率.

4.2 本文模型和其他模型对比

结果如图7所示(准确率和召回率为改进kmeans后的结果，准确率3和召回率3为使用Apriori算法的结果，横坐标单位：万条数据，纵坐标单位：百分比).

图7 本文模型和其他模型性能对比Fig.7 Comparison of the performance of this model and other models

结果分析：改进kmeans后的算法导致模型准确率和召回率提高了8%左右.对比传统关联模式挖掘的算法，本文算法具有一定的优势.说明基于转移概率的方式与挖掘时空序列轨迹相似度的方式相比，前者将一些影响轨迹趋势的别的特征也考虑进了模型中(即使没有提取相关特征)，因而会提高准确度.

5 结束语

本文基于隐马尔可夫模型，提出了一种预测移动对象下一个位置的方法—基于聚类分区域的时空轨迹预测模型，该方法用以解决时空轨迹序列的预测问题.该模型首先通过聚类将一片区域内的时空序列分成多个小区域，每个小区域内再通过聚类确定多个隐状态和发射序列，然后针对每个小区域进行隐马尔可夫模型的训练得出最终模型.预测时通过已知的时空序列，找到对应的区域模型，通过维特比算法计算出最佳隐状态序列，再结合转移矩阵做出下一个轨迹点的预测，实验表明该方法在一定程度上提高了预测准确率.

未来的工作中，可以使用机器学习或者深度学习的方式进一步提高区域的划分算法和轨迹聚类算法.