TraClus轨迹聚类算法原理及java版实现

  前一阵子我们部门接到了业务那边的一个需求。想通过用户的wifi数据计算出商场内用户最喜爱走的线路。其实说白了就是用户轨迹聚类。根据业务的需求，我们最终采用了traClus聚类算法。traClus算法相比于其它的轨迹聚类算法的一大不同点是，该算法先把一个用户的轨迹分成了若干线段，然后把基于所有用户的轨迹生成的线段放到一个集合中进行聚类。

算法本身可以划分为三个部分，分别为：

用户轨迹分段。 
分段轨迹聚类。
计算每个类簇的代表轨迹。

  虽然算法大体可以分为以上三步，但为了实现1、2两步的算法，traclus的作者提出了一种新的计算线段距离的方式。算法的第一部分是根据这种新的距离定义来对轨迹进行划分。第二部分采用了DBSCAN聚类的思想把相近的线段聚类。第3部分则找出簇中的代表轨迹（此轨迹不一般情况不会是簇中的线段直接组合而成）。下面我们详细的讲解一下距离的定义及算法的三部分。

一、线段距离的定义

  

  如上图，作者提出了三种距度量来表示两条线段间的距离。分别是垂直距离、平等距离和角度距离。需要注意的是角度距离$l_{\theta }$在$90^{\circ }\le \theta \le 180^{\circ }$。有的同学可能会说上面的公式不是对等的，线段j和线段i要怎么区分，这个我们下面会讲到。

二、traclus算法

1、用户轨迹分段

这部分的主要任务是：把用户的轨迹按照相关算法分成若干段。举个例子，有一条5个点的轨迹$p_1{p}_2{p}_3{p}_4{p}_5$,通过我们的分段算法确定p3是一个切分点，把这条轨迹分成$p_1{p}_3$ , $p_3{p}_5$两条线段。但这个切分点是怎么确定的呢，这里作者采用MDL原则（最小描述原则），通过垂直距离$d_{\perp }$与角度距离$d_{\theta }$定义了两个度量$MD{L}_{par}$和$MD{L}_{no\_par}$如果$MD{L}_{p}ar$ 大于$MD{L}_{no\_nopar}$则在当前点的前一个点进行切分。下面是具体的细节描述。

  a) 分段度量：定义了两个度量公式：
     i) $L(H)={\sum }_{j=1}^{pa{r}_i-1}lo{g}_2(len(p_{c_j},p_{c_j+1}))$
    ii) L ( D ∣ H ) = ∑ j = 1 p a r i − 1 ∑ k = c j c j + 1 − 1 { l o g 2 ( d ⊥ ( p c j p c j + 1 , p k p k + 1 ) ) + l o g 2 ( d θ ( p c j p c j + 1 , p k p k + 1 ) ) } L(D|H)=\sum^{par_i - 1} _ {j=1} \sum^{c_{j+1} - 1}_{k=c_j}\{log_2(d_{\perp}(p_{c_j}p_{c_{j+1}} , p_k p_{k+1})) + log_2(d_{\theta}(p_{c_j}p_{c_{j+1}} , p_k p_{k+1}))\} L(D∣H)=∑j=1pari​−1​∑k=cj​cj+1​−1​{log2​(d⊥​(pcj​​pcj+1​​,pk​pk+1​))+log2​(dθ​(pcj​​pcj+1​​,pk​pk+1​))}
  只看上面两个公式比较晦涩，下面我们看下在实际中是怎么应用上面两个公式的。从下图中可以看$L(H)$是直接计算两个点之间的距离然后取以2为底的对数，与上面$L(H)$的连加公式好像对不上，其实这个与我们的计算类别有关，如果我们计算的是$MD{L}_{par}$那么我们所观察的线段其实是$p_{c1}{p}_{c2}$，就是下面图中的式子，直接计算$p_{c1}{p}_{c2}($$p_1{p}_4$)的长度，然后取对数。如果我们计算的是$MD{L}_{no\_par}$，那么就是$lo{g}_2(len(p_1{p}_2)+len(p_2{p}_3)+len(p_3{p}_4))$ 。然后我们再说下$L(D|H)$这个式子，这个式子中只有垂直距离（$l_{\perp }$）和角度距离（$l_{\theta }$），我们没有看到平行距离（$l_{||}$),有同学分析说，是因为首尾相加的线段平行距离（$l_{||}$）距离为0，这个说法我不完全认同，因为从图中我们可以看出$p_1{p}_4$与$p_2{p}_3$的水平距离不是为0的。也有可能是因为加入了平行距离会使不同长度（线段数）的轨迹计算$L(D|H)$时不公平。但这只是猜测，也许作者就是不想用，或者用了发现效果不好。然后我们说另一个问题，就是在线段距离的定义的那个图中，哪条线段是$L_i$哪条是$L_j$的问题，原文中"We assign a longer line segment to $L_i$ and a shorter one to $L_j$ without losing generality." 直接简单粗暴的把较长的线段设为了$L_i$，较短的设为了$L_j$。这样做有两个好处，首先长线段为$L_i$使平行距离更容易落在$L_i$内，更好理解。而别一个好处是符合实际的物理意义，如下图，我们在算$p_4$这个切分点时，较长的线段一直是$p_1{p}_4$，只有在$p_1{p}_2{p}_3{p}_4$形成一个弧度很大的轨迹时，才会出现$p_1{p}_2$、$p_2{p}_3$、$p_3{p}_4$中某一段的距离比$p_1{p}_4$长，但一般不会出现这种情况，因为轨迹角度过大的话，路径中的某个点就成为了切分点，不会让这种大角度轨迹出现。
    
  b)切分依据：
  $MDL(cost)=L(H)+L(D|H)$
  切分的依据是上面那个公式，对于一段轨迹，我们如果要判断在某个节点是否要切分，首先要计算这个节点的$MD{L}_{par}$和$MD{L}_{no\_par}$ , 计算$MD{L}_{par}$时，我们直接按照上面的公式计算，计算$MD{L}_{no\_par}$时，只需令上面公式中的$L(D|H)$等于0，其实此时就是计算的这一段轨迹的长度。论文中给出切分依据是比较$MD{L}_{par}$和$MD{L}_{no\_par}$的大小，如果$MD{L}_{par}$ 大于 $MD{L}_{no\_par}$ 则在前一个节进行切分。 但在我们比较的是 $MD{L}_{par}$ 与 $MD{L}_{no\_par}+\xi$，其中$\xi$是一个超参。通过这个超参数，可以调节线段切分的效果。
  c)切分算法：下面是具体切分流程的伪代码。
    

2、用户轨迹聚类

这部分的主要任务是：使用DBScan的思想把第一部分划分好的线段聚成类 ，DBScan的算法不是本篇的重点，这里不做细解，过程见下图。
    

3、代表轨迹计算

  这部分的主要任务是：对上一步聚类的结果进行计算，在每一个簇中，找出一条代表轨迹。
   算法流程：
  a) 对一个簇中的所有向量（线段）求平均向量。
    
  b)把整个簇内的向量按平均向量旋转。（旋转到$x$轴平行于a)步骤求得的平均向量。
    
    
  c)使用垂直于$x$轴的sweep line延$x$轴平扫，如果与这条直线相交的向量大于等于设置的最小值（MinLns），则计算这些相交点的$y$坐标的平均值，形成点$(x_i,{\overline{y}}_i)$ , 重复此过程，直到sweep line的右边再无向量的起始或结束点。
    
    
  d)把c)步骤生成的点旋转回原来的角度，连接成一条轨迹，这条轨迹就是这个簇的代表轨迹。
  e)对所有的簇重复a)、b)、c)、d)步骤

三、算法效果图

  我实现了一版java版的traclus算法，下图使用论文中的数据(deer95)跑的一个结果，有兴趣的同学可以去下载我的源码. 跑一跑。有问题可以留言交流。