CN106383868A - 一种基于道路网络的时空轨迹聚类方法 - Google Patents

Abstract

本发明公开了一种基于道路网络的时空轨迹聚类方法，该方法包括数据采集、时空轨迹表达、时空相似性度量、子轨迹聚类和输出聚类结果步骤。通过轨迹记录设备，采集移动对象的时空轨迹数据，基于线段的轨迹表示建立时空轨迹模型，通过线性插值和语义扩充后输出轨迹文件，并对轨迹文件选取特征点，通过特征点划分子轨迹进行轨迹重构；计算子轨迹间的网络距离，作为时空相似性分析度量依据，并运用标签传播算法实现子轨迹聚类，最后输出聚类结果。本发明通过对各种时空轨迹数据进行聚类分析，可以提取时空轨迹数据中的相似性与异常特征，并有助于发现其中有意义的轨迹模式。

Description

一种基于道路网络的时空轨迹聚类方法

技术领域

本发明涉及数据挖掘技术领域，具体是一种基于道路网络的时空轨迹聚类方法。

背景技术

随着遥感、GPS、无线通讯、智能终端、物联网等技术的发展，人们可以获得大量的时空轨迹数据，通过分析恐怖分子的活动路径来更好的指导反恐活动；通过研究城市中汽车的移动轨迹，对城市道路规划做出改善，以缓解城市拥堵问题；通过对飓风的移动轨迹数据进行研究，找到其移动模式，以便在賜风来临之前做好防灾工作，动物的迁徙轨迹，智能交通管理，生物学领域中精子运动的轨迹和人的笔迹轨迹，化学领域中分子的运动轨迹，体育领域中球员的运动轨迹等等。越来越多的轨迹数据被应用于时间地理学、交通地理学、生态地理学、医学等相关研究，这些轨迹中蕴含着大量有价值的信息与知识。

智能交通系统建立了巨大的数据库资源，在这些海量数据中隐含地存在着大量对路径规划有用的知识。例如在交通导航系统中，根据对轨迹的相似性查询技术可以得到历史相似运动轨迹，基于这些得到的知识就可以推荐一组比较合适的行车路线或估计所需要的时间给新手或迷路司机。除此之外，通过轨迹的相似性查询还可以发现频繁发生交通堵塞的区域中大部分车辆的运动模式，根据这些结果可以规划开发一些新的道路或者为出行的车辆提供路况信息等建议。聚类分析作为数据挖掘技术的一项任务，是一种数据预处理过程，是进一步分析和处理数据的基础。基于道路网络空间的移动对象轨迹聚类技术可为路径规划的合理性与最优性提供有力的保障，为用户提供相应的决策支持，轨迹聚类的研究在实际应用中具有重要意义。

移动对象的轨迹聚类从其产生至今，就一直受到各界研究人员的关注。采用轨迹聚类方法，通过发现相似的运动轨迹、提取运动特征模式等方式，可以发现移动对象的运动规律和行为模式。时空轨迹是移动对象的位置和时间的记录序列。时空轨迹数据不同于聚类分析中传统数据的数据结构，聚类分析在度量传统数据之间的相似性时，大多考虑的是数据点和数据点之间的相似性，而轨迹是由点的序列构成的，不同轨迹之间的点有的离得非常近有的离得非常远。在考虑轨迹之间的相似性时就不能只考虑点和点之间的相似性了，而要考虑由点构成的线与线之间的相似性。

移动对象的时空轨迹聚类算法研究，主要方法目前可分为两类：一类是以整条轨迹为分类对象进行分类，优点在于其在评价轨迹对象之间的相似性方面比较直观，用户输入参数的影响比较小；但缺点是对于复杂的轨迹，此类方法无法分辨出轨迹的局部特征，且对高维度的轨迹数据的聚类效果不好。另一类是将复杂轨迹按照某种标准分割为相对平滑的轨迹段，然后以分割得到的轨迹段为分类对象进行分类，优点在于可以很好地识别轨迹的局部特征，有效地处理高维度的轨迹数据，结合基于密度的聚类方法，能够很好地发现任意形状的轨迹段聚类，但缺点是轨迹分割方法对聚类结果有较大的影响，不少专家学者提出了许多不同的分割方法，适用于不同的应用领域，但领域知识对聚类结果的影响较大。

尽管如此，时空轨迹聚类方法仍面临许多困难与挑战，有待进一步研究和解决，主要体现在以下5个方面：

(1)当前大部分时空轨迹数据聚类方法仍然是将时间看作原空间对象的附加维，这种处理方式难免使时间与空间有所分隔，与人们对于事物的直观认识有出入；

(2)现有聚类方法对于某些轨迹数据类型并不完全适用，例如人口迁移轨迹数据，该类数据的时间维是不等长的，但是由于具有年龄的语义，时间维不能拉伸，这种情况目前还没有很好的方法可以处理；

(3)聚类结果在转换成知识的过程中存在一些问题，例如所发现的知识或者过于简单，近乎于常识，或者过于复杂，让人们无法直观理解；

(4)海量的轨迹数据的不断产生，一方面为研究者提供了丰富的数据源，但另一方面，也要求研究者从中选择有效的数据并提高算法效率；

(5)在处理与人有关的时空轨迹时，如何保护对象的隐私等也成为研究者应当考虑的问题。

本发明采用标签传播算法进行轨迹聚类，并结合道路网络应用领域，在局部算法上进行突破和改进，并从总体上提出了一种基于道路网络的时空轨迹聚方法。

目前大多数空间聚类算法，主要针对几何空间数据对象，而在交通领域，空间对象的访问受限于道路网络，很多兴趣点之间的路网距离和欧式距离差别很大。因此，在交通领域，研究基于路网拓扑的聚类分析算法，更具有应用价值。

常见的聚类算法采用欧氏距离对数据对象的相似度进行度量，通常发现的是一些密度和大小比较接近的球状的或圆状的聚类。但在现实世界中聚类不只是球状或圆状，可以是任意形状，因此设计能够发现任意形状聚类的算法具有非常重要的现实意义。标签传播算法(LPA)可以通过相近节点之间的标签的传递来学习分类，它不受数据分布形状的局限，可以克服一些算法只能发现“类圆形”结构的缺点。只要同一类的数据在空间分布上是相近的，那么不管数据分布是什么形状，都能通过标签传播将它们分到同一个类里。

同时由于hansdorff distance依赖于参数ε，而且在使用欧几里得距离删减流量簇和计算hausdorff邻域这两个过程中都使用到了ε，所以最终聚类结果的好坏与参数的选择密切相关，然而参数的确定往往依靠经验或专家知识，显然不利于用户的使用，同时也使得聚类分析的结果难以得到有效的保证。

聚类的伸缩性要求主要是指算法不仅能够处理较小的数据集合，还要能够处理含有百万、千万条数据记录的海量数据对象。为了便于扩展，就要求算法的执行效率不能过低，最好能够在线性时间或者多项式时间内完成。当轨迹数据增多，使用hausdorffdistance计算量就会增大，效率迅速降低。然而标签传播算法(LPA)只需利用少量的训练标签指导，利用未标注数据的内在结构、分布规律和邻近数据的标记，即可预测和传播未标记数据的标签，然后合并到标记的数据集中。该算法操作简单、运算量小，适合大规模数据信息的挖掘和处理，准确性高于基于hausdorff distance的聚类方法。

发明内容

本发明的目的在于克服现有技术的不足，基于经典聚类算法进行改进，提供一种基于道路网络的时空轨迹聚类方法。

本发明的目的是通过以下技术方案来实现的：一种基于道路网络的时空轨迹聚类方法，它包括以下步骤：

S1：数据采集，利用轨迹记录设备采集移动对象的时空轨迹数据，或者以用户真实采集的GPS数据为基础，作为数据输入；

S2：时空轨迹表达，基于局部插值模型，通过线性插值和语义扩充建立模型输出轨迹文件，并对轨迹文件选取特征点，通过特征点，基于速度和角度划分子轨迹进行轨迹重构；

S3：时空相似性度量，计算子轨迹间的网络距离，作为时空相似性度量依据；

S4：子轨迹聚类，运用标签传播算法，将已标记的数据作为标签传播的源，并保持一致的标签，准确实现子轨迹聚类；

S5：输出聚类结果。

所述的时空轨迹表达S2包括如下子步骤：

S21：基于线段的轨迹表示建立时空轨迹模型，在步骤S1输入的时空轨迹数据中，提取停留点进行语义扩充，并将停留点轨迹转变为位置轨迹；

S22：基于局部插值模型的时空轨迹数据表达，将时空轨迹网格离散化后，利用反距离权重法计算落在每个网格内轨迹段的各个属性的值，根据临近的采样点属性对轨迹段进行插值；

S23：基于速度和角度轨迹段分离法寻找特征点，并加入速度约束来分割轨迹；

S24：利用基于网路的子轨迹划分算法，输入完整轨迹和道路交汇点，输出重新划分子轨迹后的完整轨迹和划分特征点集合。

所述的寻找特征点S23包括如下子步骤：

S231：在给定距离内，首先扫描计算每个点，判断出候选关键点和候选变速点；

S232：清除多余候选关键点和多余候选变速点，判断出关键点和变速点，利用关键点确定轨迹段；

S233：首先，开始点被加入到特征点集合；

S234：进入循环，找出轨迹特性快速改变的点并加入到特征点集合中，分别用变量表示出相邻子段的方向和速度变化，并分别用变量表示出累计的方向和速度变化；

S235：判断累计变量，若超过阈值，则前一个点被选为特征点。

所述的子轨迹划分算法S24包括如下子步骤：

S241：首先输入完整轨迹和道路交汇点，将移动对象的第一个采样点加入到特征点集合中；

S242：接下来循环判断每个子轨迹的时空起点和时空终点是否被道路交汇点隔开；

S243：如果隔开则重新划分当前子轨迹并加入整轨迹，然后重置当前子轨迹，继续循环计算该子轨迹的特征点；

S244：当子轨迹段的时空起点和时空终点在一条道路段中时，判断子轨迹段的方向、速度和地理空间变化，找出轨迹特性快速改变的点并加入到特征点集合中；

S245：如果该点不能作为特征点，则将当前子轨迹合并，并将合并后的子轨迹重新加入特征点的计算中；

S246：输出重新划分子轨迹后的完整轨迹和划分特征点集合。

所述的时空相似性度量S3包括如下子步骤：

S31：判断特征点A、B是否在同一路段；

S32：若特征点A、B在同一路段，则利用A、B两点在二维平面的坐标计算两点之间的距离；

S33：若特征点A、B不在同一路段，则寻找A、B两点的相邻交叉口n_a、n_b，分别计算A、B两点到各自相邻交叉口的距离，再计算相邻交叉口之间的最短距离，求和即为特征点A、B之间的网络距离。

步骤S33所述的相邻交叉口之间最短距离d_ij的计算，包括如下子步骤：

S331：先扫描与n_a相邻的交叉路口并同时记录此路段的长度；

S332：在扫描过程中记下每条路径所经过的总长度，选取距离值最小的一条；

S333：在遍历的过程中，判断循环提前结束的条件；

S334：若已找出一条到达n_b的路径，但是当前遍历没有结束，且到目前为止的总代价已经大于先前检索出的路径，则此次循环提前结束，进入下一条可能路径的查询；

S335：否则继续执行遍历过程，每条相邻路段都遍历一次，直到到达路口n_b，即找到相邻交叉口之间的最短路径。

所述的子轨迹聚类S4包括如下子步骤：

S41：计算边的权重w_ij，将所有子轨迹段作为节点，创建一个完全连接图，计算公式如下：

$w_{ij}=\exp \left(\frac{{d_{ij}}^2}{{\mathrm{\σ}}^2}\right)=\exp \left(-\frac{\underset{d}{\overset{D}{\Σ}}{(x_i^d-x_j^d)}^2}{{\mathrm{\σ}}^2}\right)$

其中，d_ij表示两个子轨迹段之间的网络距离，即相应特征点之间最短距离，σ为自行定义的参数，x_i表示第i节点位置的横坐标，x_j表示第j节点位置的横坐标；

S42：定义参数σ，基于kruskal算法找到树中第一个连接不同标签的节点的边，以这条边的长度d₀作为不同类之间最小距离的一个启发值，设置σ＝d₀/3，因此这条边的权重为0，从而保证在类内标签传播阻碍为0，以子轨迹段之间的最短距离作为不同类之间最小距离的一个启发值；

定义一个(l+u)×(l+u)的概率传递矩阵T：

$T_{ij}=P(j\→i)=\frac{w_{ij}}{\underset{k=1}{\overset{l+u}{\Σ}}{w}_{kj}}$

其中，T_ij表示节点j到i的传播概率；

定义一个(l+u)×C的标注矩阵Y，令Y_ic＝δ(y_i,c)，它的第i行代表着节点y_i的标注概率，第c列代表类别，若Y_ic＝1，则表示节点y_i是属于c类别，若Y_ic＝0则通过概率传递，使其概率分布集中于给定类别，然后通过边的权重值来传递节点标签，保证标注矩阵Y的每行都是标准化的；

S43：传播标签Y←TY，所有节点将其标签传播到所有的其他节点；

S44：对Y进行归一化，保持类标签的概率表达；

S45：重置已标记数据的类标签，重复步骤S43，将已标记的数据能够作未标签传播的源，保持一致的标签，直到矩阵Y收敛。

本发明的有益效果是：运用本设计方法，降低了轨迹的存储空间，提高了轨迹的时空相似性搜索效率和准确性，能够有效发现时空邻近的轨迹，且具有较好的可伸缩性，具有极大的应用前景。

附图说明

图1为本方法数据流向示意图；

图2为本方法开放角示意图；

图3为本方法变速点示意图；

图4为本方法候选关键点的计算示意图；

图5为本方法多余候选关键点的清除示意图；

图6为本方法局部插值模型的示意图；

图7为本方法时空路径可视化方法示意图；

图8为本方法寻找特征点的流程图；

图9为本方法子轨迹划分的流程图；

图10为本方法分段线性插值的示意图；

图11为本方法计算网络距离的示意图；

图12为本方法基于标签传播算法的子轨迹聚类流程图。

具体实施方式

下面结合附图进一步详细描述本发明的技术方案，但本发明的保护范围不局限于以下所述。

如图1为本发明的系统框图，首先建立时空轨迹模型：

定义1坐标点：一个带时间戳的坐标点定义为P＝(Lat,Lngt,T)，其中Lat为纬度，Lngt为经度，T为时间戳。

定义2GPS轨迹：一条GPS轨迹Tra定义为按时序组成的一组坐标点序列，Tra＝P1->P2->…～pn，其中pi.T<pi+1.T。

定义3停留点川：给定时间阈值θ_t和距离阈值θ_d，停留点s定义为1个时间段内用户轨迹中一段连续坐标点P＝(P_m,P_m+1,…,P_n)，其中：

$\&ForAll;m<i\&le;n,Dis\tan ce(p_m,p_i)\&le;{\mathrm{\&theta;}}_d,|p_n,T-p_m.T|\&GreaterEqual;{\mathrm{\&theta;}}_t$

将停留点s的相关属性记为(s.Loat,s.Lngt,s.arvT,s.levT)，其中s.Loat和s.Lngt表示P中各坐标经度和纬度的平均值，s.arvT＝pm.T为用户到达停留点的时间，s.levT＝pn.T为离开停留点的时间。相对于坐标点，停留点含有更重要的语义信息，可能代表餐馆、商店、电影院、旅游景点等。

定义4停留点轨迹：一段按时序组织的停留点序列其中Δt₁＝s_i+1.arvT-S_i.levT表示停留点之间的行程时间。

定义5位置轨迹：若给定一停留点轨迹为位置轨迹定义为其中s_i∈c对所有停留点作聚类分析后，每个停留点s_i都有所属簇的标号c_i，把簇的序列称为用户的位置轨迹。

以用户真实采集的GPS轨迹数据为基础，定义相关变量，提取出停留点，进行语义扩充，并通过对停留点进行聚类分析生成基础聚类，把停留点轨迹转变为位置轨迹。

然后对时空轨迹进行网格离散化：假设相邻记录点间对象作匀速直线运动，该规则可以用线性插值方法表达，局部插值模型如图6，采用时空路径(Space-time Path)的方式来可视化表达如图7，实线表示的是时空路径，虚线为时空路径在空间维上的投影。这种表达方式将二维的空间和一维的时间整合到一个三维坐标系中表示，每个记录点的x、y坐标对应记录点的空间坐标，第三维坐标则对应记录点的时间值。

时空轨迹(Trajectory，简称Tr)定义为由若干个多维度的采样点按照时间顺序排列组成的一个序列，Tr_i＝{P₁，P₂，…，Pm}(1≤i≤m)，其中P为一个多属性采样点，表示为：Pj＝{Loc，T，A₁，…，A_i，…，A_m}。其中，A_i表示采样点P_j的第i个属性，Loc代表采样点Pj的空间位置，是一个多维度的位置点，T表示采样点P_j的采样时间。移动对象在轨迹中除采样点外任意一点的属性信息，可以通过一定的规则，根据其所在轨迹段的前后邻近采样点的属性做近似的插值得到。

A_ji＝fi(A_(m-1)i，…，A_mi，A_ni，…，A_(n+h)i，dis_(m-1)i，…，dis_mi，dis_ni，…，dis_(n+h)i)。其中，A_ji表示P_j点的第i个属性，A_mi表示P_j点之前相邻采样点的属性，A_ni表示P_j点之后相邻采样点的属性，dis_mi表示P_j点与之前相邻采样点的距离，dis_ni表示P_j点与之后相邻采样点的距离，并且l与h之和等于所设定的用来插值的采样点数量s。fi表示根据P_j前后采样点属性Ai和P_j与其他采样点之间距离，对P_j的Ai属性进行插值的方法。

轨迹被离散后，利用反距离权重法计算落在每个网格内轨迹段的各个属性的值，根据临近的采样点属性对轨迹段进行插值。用户可以设置距离在插值过程中对属性变化的影响大小，反距离权重法通用方程如下：

$Z_o=\underset{i=1}{\overset{s}{\&Sigma;}}{Z}_i\frac{1}{d_i^k}/\underset{i=1}{\overset{s}{\&Sigma;}}\frac{1}{d_i^k}$

式中：Z_o为点O的估计值；Z_i为同一轨迹其他已知采样点的Z值；d_i为采样点i与点O之间的空间距离；s为在估算中用到的控制点的数目；k为指定的幂。

基于速度和角度划分轨迹段，如图2所示，在给定距离(dmin,dmax)内，如果(P-,P,P+)三点构成的夹角小于或等于给定值αmax，该角称为开放角，则就可以判定该P点为关键点；如果在P点邻域距离(dmin,dmax)内，任意点与P点的速度差大于或等于ε₂，则不管P点的开放角是多大，都判定P为关键点，该点称为变速点，如图3所示，两个相邻关键点间轨迹点的连线集合就构成了轨迹段。

轨迹上的点如果曲度很大，那么它形状突出，则它在这条轨迹上的位置比较关键，本方法认为点速度变大的在轨迹上也扮演重要角色。通过开放角和变速点来确定关键点能比较好的对轨迹分段，其中开放角考虑到轨迹段的距离和方向变化，变速点只考虑速度的变化，加入速度约束来分割轨迹，关键点的检测具体如图4，图5所示，该算法分2步：

第1步为扫描计算每个点，如图4所示，如果满足：d_min ²≤|p-p^-|²≤d_max ²，d_min ²≤|p-p⁺|²≤d_max ²，α≤α_max这3个条件，那么该点就为候选开放角关键点。或者如果它们的速度满足：那么该点是候选变速关键点。

第2步为清除多余候选关键点，如图5所示，如果α(p)＞α(p_v)，并且pv在p的邻域内，即|p-p_v|²≤d_max ²(或者也可以给条件为|p-p_v|²≤d_min ²)，那么就删除p点，保留pv点；但是如果该p点同时已经被判定为候选变速点，那么保留p点，删除pv点；而对于候选速度关键点，如果该点同时已经被判定为候选开放角，那么就保留该点；如果pv、p都为候选变速点，而且pv在p的邻域内，若α(p)＞α(p_v)，则删除p点，保留pv点。

如图8所示，进行轨迹分割的目的是要找出轨迹的速度和方向显著变化的点，称之为特征点。轨迹分割成子段的算法T-part描述如下：

输入：(1)轨迹TR＝{<pi，ti>(|1≤i≤n)}

(2)4个阈值参数：Min_dc1，Min_vc1，Min_dc2，Min_vc2；

输出：特征点的集合CP。

01：Add<p1，t1>into the set CP；

02：StarIndex＝1；length＝1；

03：while StartIndex+length≤n do{

04：curr＝StarIndex+length；

05：compute direction changes dc1 and speed changes vc1 between linesegments pcurr-2pcurr-1 and pcurrI-1pcurr

06：compute direction changes dc2 and speed changes vc2 between linesegments pstartIndexpstartIndex+1 and pstartIndexpcurr

07：if(dc1>Min_dc1)or(sc1>Min_sc1)or(dc2>Min_dc2)or(sc2>Min_sc2)

08：{Add Pcurr-1 into the set CP

09：StartIndex＝curr-1；length＝1；}

10：else

11：length＝length+1；}

12：Add pninto the set CP；

首先，开始点p1被加入到特征点集合CP中(第1行)，接下来的循环中找出轨迹特性快速改变的点加入到CP中(第3～11行)，dc1和vc1分别表示相邻子段的方向和速度变化，dc2和vc2分别表示累计的方向和速度变化，尽管有时相邻变化不大，累计的变化却不小。如果变化量超过阈值，则前一个点pcurr-1被选为特征点，算法的时间复杂度是O(logn)。

针对城市交通网网路轨迹的划分，子轨迹划分算法流程示意图，如图9所示，将整轨迹划分程子轨迹的算法描述如下：

算法：子轨迹划分算法SubTR_Partition

输入：完整轨迹TR＝{idTR,subTr1,…,subTrn}和道路交汇点Crossing；

输出：重新划分子轨迹后的TR和划分特征点集合CPSet。

具体算法过程：

首先，将移动对象的第一个采样点P1加入到特征点集合CPSet中，接下来的循环判断每个子轨迹的时空起点和时空终点是否被道路交汇点Crossing隔开，如果隔开则将当前子轨迹subTri{P1,P2}重新划分成subTri{P1,Crossing}加入整轨迹tmpTR，然后重置当前子轨迹为subTri{Crossing,P2}(6-11行)，继续循环计算该子轨迹的特征点；当子轨迹段的时空起点和时空终点在一条道路段中时，判断子轨迹段的方向、速度和地理空间变化，将找出轨迹特性快速改变的点加入到CP中，如果该点不能作为特征点，则将subTri-1和subTri合并，并将合并后的子轨迹重新加入特征点的计算。

所述时空相似性度量，首先计算任意子轨迹段之间的网络距离，即求相应特征点之间的网络距离。例如计算特征点A、特征点B之间的网络距离，将A、B位置信息用坐标分别表示为(x_A,y_A)，(x_B,y_B)表示，它们之间的网络距离分为两种情况：

1、如果A、B在同一路段上，Rid_A＝Rid_B时，

2、如果A、B不在同一路段，且n_a和n_b分别是A和B相邻的交叉口，则dist(A,B)＝dist(A,n_a)+dist(n_a,n_b)+dist(n_b,B)，其中dist(A,n_a)和dist(n_b,B)都是在求同一路段上两点之间的距离，dist(n_a,n_b)是交叉口n_a到n_b的最短距离。求n_a到n_b最短距离:先扫描与n_a相邻的交叉路口并同时记录此路段的长度，然后向外延伸，每条相邻路段都遍历一次，直到到达路口n_b。在这个过程中记下每条路径所经过的总长度，选取值最小的一条。在遍历的过程中，会有一个循环提前结束的条件，就是如果已找出一条到达n_b的路径，但是当前遍历没有结束，且到目前为止的总代价已经大于先前检索出的路径，则此次循环提前结束，进入下一条可能路径的查询。这种方法不需要遍历整个路网，也不需要知道所有边的权值，减少了计算时间。

例：如图11所示，计算P3到P7路网距离，兴趣点P7不是交叉点，则找到与其相邻的交叉点P6和P9，并记录P7到这两点的距离，然后从P3出发，分别查找到P6和P9的最短路径。首先扫描P3的三个相邻交叉点P2，P5，P10，并记录三段的距离，然后对这三个点循环执行扫描相邻交叉点的操作，直到找到P6或P9为止。在这个过程中完成了P3-P5-P6的查找，并发现其他路径还没有找到P6或P9并且长度已经大于该路径的长度，则停止其他路径的扫描，得出dist(p3,p7)＝dist(p3,p5)+dist(p5,p6)+dist(p6,p7)。

作为较佳的实施例之一，这里提供有一种简化方法：在遍历路网时，不是遍历所有的相邻边，而是线判断点B在点A的什么方向(通过坐标值判断)，上，下，左，右，左上，左下，右上，右下，确定相对方向后就只走相应方向的边。例如若B在A的右下方，则从n_a出发，只向右和向下走，其他相邻边则忽略，这样能减少遍历及计算的次数，而且在大部分情况下，这种方法查询的结果是正确的。

运用标签传播算法实现子轨迹聚类，如图12所示，概念定义如下：

令(x₁,y₁)…(x_l,y_l)是已标注数据，Y_L＝{y₁…y_l}∈{1…C}是类别标签，类别数C已知，且均存在与标签数据中。令(x_l+1,y_l+1)…(x_l+u,y_l+u)为未标注数据，Y_U＝{y_l+1…y_l+u}不可观测，l＜＜u，令数据集X＝{x₁…x_l+u}∈R^D。问题转换为：从数据集X中，利用Y_L的学习，为未标注数据集Y_U的每个数据找到对应的标签。首先计算任意两个特征点之间的网络距离，详细步骤：

S31：判断特征点A、B是否在同一路段；

S32：若特征点A、B在同一路段，则利用A、B两点在二维平面的坐标计算两点之间的距离；

S33：若特征点A、B不在同一路段，则寻找A、B两点的相邻交叉口n_a、n_b，分别计算A、B两点到各自相邻交叉口的距离，再计算相邻交叉口之间的最短距离，求和即为特征点A、B之间的网络距离。

步骤S33所述的相邻交叉口之间最短距离d_ij计算，包括如下子步骤：

S331：先扫描与n_a相邻的交叉路口并同时记录此路段的长度；

S332：在扫描过程中记下每条路径所经过的总长度，选取距离值最小的一条；

S333：在遍历的过程中，判断循环提前结束的条件；

S334：若已找出一条到达n_b的路径，但是当前遍历没有结束，且到目前为止的总代价已经大于先前检索出的路径，则此次循环提前结束，进入下一条可能路径的查询；

S335：否则继续执行遍历过程，每条相邻路段都遍历一次，直到到达路口n_b，即找到相邻交叉口之间的最短路径。

计算出子轨迹段之间的网络距离以后，再运用标签传播算法进行子轨迹聚类，详细步骤：

S41：计算边的权重w_ij，将所有子轨迹段作为节点，创建一个完全连接图，计算公式如下：

$w_{ij}=\exp \left(\frac{{d_{ij}}^2}{{\mathrm{\&sigma;}}^2}\right)=\exp \left(-\frac{\underset{d}{\overset{D}{\&Sigma;}}{(x_i^d-x_j^d)}^2}{{\mathrm{\&sigma;}}^2}\right)$

其中，d_ij表示两个子轨迹段之间的网络距离，即相应特征点之间最短距离，σ为自行定义的参数，x_i表示第i节点位置的横坐标，x_j表示第j节点位置的横坐标；

S42：定义参数σ，基于kruskal算法找到树中第一个连接不同标签的节点的边，以这条边的长度d₀作为不同类之间最小距离的一个启发值，设置σ＝d₀/3，因此这条边的权重为0，保证在类内标签传播阻碍为0，以子轨迹段之间的最短距离作为不同类之间最小距离的一个启发值；

定义一个(l+u)×(l+u)的概率传递矩阵T：

$T_{ij}=P(j\&RightArrow;i)=\frac{w_{ij}}{\underset{k=1}{\overset{l+u}{\&Sigma;}}{w}_{kj}}$

其中，T_ij表示节点j到i的传播概率；

定义一个(l+u)×C的标注矩阵Y，令Y_ic＝δ(y_i,c)，它的第i行代表着节点y_i的标注概率，第c列代表类别，若Y_ic＝1，则表示节点y_i是属于c类别，若为Y_ic＝0则通过概率传递，使其概率分布集中于给定类别，然后通过边的权重值来传递节点标签，保证标注矩阵Y的每行都是标准化的；

S43：传播标签Y←TY，所有节点将其标签传播到所有的其他节点；

S44：对Y进行归一化，保持类标签的概率表达；

S45：重置已标记数据的类标签，重复步骤S43，将已标记的数据能够作未标签传播的源，保持一致的标签，直到矩阵Y收敛。

程序结束，输出聚类结果。

Claims (7)

1.一种基于道路网络的时空轨迹聚类方法，其特征在于，它包括以下步骤：

S1：数据采集，利用轨迹记录设备采集移动对象的时空轨迹数据，或者以用户真实采集的GPS数据为基础，作为数据输入；

S2：时空轨迹表达，基于局部插值模型，通过线性插值和语义扩充建立模型输出轨迹文件，并对轨迹文件选取特征点，通过特征点，基于速度和角度划分子轨迹进行轨迹重构；

S3：时空相似性度量，计算子轨迹间的网络距离，作为时空相似性度量依据；

S4：子轨迹聚类，运用标签传播算法，将已标记的数据作为标签传播的源，并保持一致的标签，准确实现子轨迹聚类；

S5：输出聚类结果。

2.根据权利要求1所述的一种基于道路网络的时空轨迹聚类方法，其特征在于，所述的时空轨迹表达包括如下子步骤：

S21：基于线段的轨迹表示建立时空轨迹模型，在步骤S1输入的时空轨迹数据中，提取停留点进行语义扩充，并将停留点轨迹转变为位置轨迹；

S22：基于局部插值模型的时空轨迹数据表达，将时空轨迹网格离散化后，利用反距离权重法计算落在每个网格内轨迹段的各个属性的值，根据临近的采样点属性对轨迹段进行插值；

S23：基于速度和角度轨迹段分离法寻找特征点，并加入速度约束来分割轨迹；

S24：利用基于网路的子轨迹划分算法，输入完整轨迹和道路交汇点，输出重新划分子轨迹后的完整轨迹和划分特征点集合。

3.根据权利要求2所述的一种基于道路网络的时空轨迹聚类方法，其特征在于，所述的寻找特征点包括如下子步骤：

S231：在给定距离内，首先扫描计算每个点，判断出候选关键点和候选变速点；

S232：清除多余候选关键点和多余候选变速点，判断出关键点和变速点，利用关键点确定轨迹段；

S233：首先，开始点被加入到特征点集合；

S234：进入循环，找出轨迹特性快速改变的点并加入到特征点集合中，分别用变量表示出相邻子段的方向和速度变化，并分别用变量表示出累计的方向和速度变化；

S235：判断累计变量，若超过阈值，则前一个点被选为特征点。

4.根据权利要求2所述的一种基于道路网络的时空轨迹聚类方法，其特征在于，所述的子轨迹划分算法包括如下子步骤：

S241：首先输入完整轨迹和道路交汇点，将移动对象的第一个采样点加入到特征点集合中；

S242：接下来循环判断每个子轨迹的时空起点和时空终点是否被道路交汇点隔开；

S243：如果隔开则重新划分当前子轨迹并加入整轨迹，然后重置当前子轨迹，继续循环计算该子轨迹的特征点；

S244：当子轨迹段的时空起点和时空终点在一条道路段中时，判断子轨迹段的方向、速度和地理空间变化，找出轨迹特性快速改变的点并加入到特征点集合中；

S245：如果该点不能作为特征点，则将当前子轨迹合并，并将合并后的子轨迹重新加入特征点的计算中；

S246：输出重新划分子轨迹后的完整轨迹和划分特征点集合。

5.根据权利要求1所述的一种基于道路网络的时空轨迹聚类方法，其特征在于，所述的时空相似性度量包括如下子步骤：

S31：判断特征点A、B是否在同一路段；

S32：若特征点A、B在同一路段，则利用A、B两点在二维平面的坐标计算两点之间的距离；

S33：若特征点A、B不在同一路段，则寻找A、B两点的相邻交叉口n_a、n_b，分别计算A、B两点到各自相邻交叉口的距离，再计算相邻交叉口之间的最短距离，求和即为特征点A、B之间的网络距离。

6.根据权利要求5所述的一种基于道路网络的时空轨迹聚类方法，其特征在于，步骤S33所述的相邻交叉口之间最短距离d_ij的计算，包括如下子步骤：

S331：先扫描与n_a相邻的交叉路口并同时记录此路段的长度；

S332：在扫描过程中记下每条路径所经过的总长度，选取距离值最小的一条；

S333：在遍历的过程中，判断循环提前结束的条件；

S334：若已找出一条到达n_b的路径，但是当前遍历没有结束，且到目前为止的总代价已经大于先前检索出的路径，则此次循环提前结束，进入下一条可能路径的查询；

S335：否则继续执行遍历过程，每条相邻路段都遍历一次，直到到达路口n_b，即找到相邻交叉口之间的最短路径。

7.根据权利要求1所述的一种基于道路网络的时空轨迹聚类方法，其特征在于，所述的子轨迹聚类包括如下步骤：

S41：计算边的权重w_ij，将所有子轨迹段作为节点，创建一个完全连接图，计算公式如下：

$w_{ij}=\exp \left(\frac{{d_{ij}}^2}{s^2}\right)=\exp (-\frac{\underset{d}{\overset{D}{\&Sigma;}}{(x_i^d-x_j^d)}^2}{s^2})$

其中，d_ij表示两个子轨迹段之间的网络距离，即相应特征点之间最短距离，σ为自行定义的参数，x_i表示第i节点位置的横坐标，x_i表示第j节点位置的横坐标；

S42：定义参数σ，基于kruskal算法找到树中第一个连接不同标签的节点的边，以这条边的长度d₀作为不同类之间最小距离的一个启发值，设置σ＝d₀/3，因此这条边的权重为0，从而保证在类内标签传播阻碍为0，以子轨迹段之间的最短距离作为不同类之间最小距离的一个启发值；

定义一个(l+u)×(l+u)的概率传递矩阵T：

$T_{ij}=P(j\&RightArrow;i)=\frac{w_{ij}}{\underset{k=1}{\overset{l+u}{\&Sigma;}}{w}_{kj}}$

其中，T_ij表示节点j到i的传播概率；

定义一个(l+u)×C的标注矩阵Y，令Y_ic＝δ(y_i,c)，它的第i行代表着节点y_i的标注概率，第c列代表类别，若Y_ic＝1，则表示节点y_i是属于c类别；若Y_ic＝0则通过概率传递，使其概率分布集中于给定类别，然后通过边的权重值来传递节点标签，保证标注矩阵Y的每行都是标准化的；

S43：传播标签Y←TY，所有节点将其标签传播到所有的其他节点；

S44：对Y进行归一化，保持类标签的概率表达；

S45：重置已标记数据的类标签，重复步骤S43，将已标记的数据能够作未标签传播的源，保持一致的标签，直到矩阵Y收敛。