0 引言

近年来随着移动定位服务和云计算技术的发展，搜集和处理移动定位信息已成为现实。由移动定位数据组成的海量数据库推动了关于移动定位数据研究的发展。作为典型的移动定位数据之一，时空轨迹数据是指沿时间顺序对空间移动对象的数值型记录，包括采样点位置、采样时间、速度等，而空间对象的这些特征都有可能随着时间的推移而发生变化，人们尝试着将这些可感而不可知的感性知识转化为有现实意义的数字信息，用于生活，提供便利。然而面对规模庞大的时空轨迹数据，如何快速地从中挖掘有用价值是人们必须面对的一个重要问题。

为了从时空轨迹数据中获取有用信息，聚类技术被用于对空间移动对象时空特征的提取。目前研究重点主要集中在轨迹数据模型定义、相似性度量和轨迹聚类方法三个方面。文献[1]基于全区间模型对轨迹进行聚类，提出一种基于采样和密度的轨迹聚类算法，使用轨迹全体采样点来度量整条轨迹的相似性，并采用基于密度的聚类算法完成聚类。文献[2]根据已有时空轨迹数据样本提取轨迹线索特征并进行汇总，以线索相似性来衡量两条轨迹之间的距离。然而将轨迹作为整体进行聚类的方法可能会忽略存在于整体相似度低的轨迹之间的局部相似性。文献[3]提出了一种划分-归组(Partition-and-Group)聚类框架，该方法首先使用最小描述长度 (Minimum Description Length, MDL)对轨迹进行子轨迹段划分，然后使用基于密度的聚类方法(Density-Based Spatial Clustering of Applications with Noise, DBSCAN)对划分后的子轨迹段进行聚类。其中MDL方法产生于对某种模型进行编码压缩的要求，其原理是计算出总描述长度最小的模型，应用到轨迹划分中时，可能会导致轨迹划分不具有实际代表性，失去其现实意义。文献[4]首先根据相邻轨迹样本点的转向角对轨迹进行划分，其次通过计算轨迹段的结构(方向、速度、转角、位置)相似度来判断轨迹的匹配程度，进而完成轨迹聚类。该提取特征点的方法适用于相邻样本点角度变化较明显的如道路轨迹数据中, 然而在相邻样本点持续局部小角度变化过程中形成整体大角度变化的情况下，该判断特征点的方法可能会导致移动轨迹的划分失去其准确性。

时空轨迹数由位置坐标和相应的时间标记合成，因此，时间特征对分析移动物体的运行模式非常关键。然而现有轨迹聚类算法中常常忽略了对时间特征的相似性度量，如文献[1-4]，为此，本文提出基于时空模式的轨迹数据聚类方法，以划分-归组聚类框架为基础，从轨迹特征点提取、时空相似性度量和轨迹聚类算法三个方面进行改进。首先为解决MDL或“相邻样本点转角”方法提取轨迹特征点有所不足的问题，借鉴图像曲线边缘信息检测的方法提取轨迹特征点，生成具有代表性的子轨迹段集合；其次采用DBSCAN算法同时根据空间和时间相似性对子轨迹段进行聚类，并使用搜索网格提高计算效率。

1 轨迹特征点的提取 1.1 相关定义

时空轨迹数据记录了持续移动对象在某时刻出现在相应位置上的一系列信息，是其运动状态下时空属性的表现。参考文献[3]对时空轨迹数据的定义进行扩展，在这里给出时空轨迹数据集、特征点集合以及子轨迹集合的数学定义。

定义1 时空轨迹数据集。给定时空轨迹数据集TRS，共有tranum条轨迹，表示为TRS={TR₁, TR₂, …, TR_tranum}，集合中的每条轨迹都由若干个多维样本点组成，任意一条轨迹TR_i表示为：TR_i={P₁, P₂, …，P_itnum}，其中1≤i≤tranum，itnum为第i条时空轨迹样本点总数。轨迹中任意一个样本点P_j表示为P_j=〈p_num, p_traid, p_x, p_y, p_time〉，显示出在p_time时刻，轨迹p_traid中样本点p_num的位置为(p_x, p_y)，其中1≤j≤itnum。

定义2 特征点集合。对第i条轨迹进行子轨迹段划分，得到能够代表该条轨迹的特征点集合CHA_TR_i，表示为CHA_TR_i={tc₁, tc₂, …, tc_icnum}，其中icnum为第i条轨迹特征点总数，集合中tc_i为该条轨迹的特征点，tc_i定义如样本点。

定义3 子轨迹段集合。根据所有轨迹的特征点集合生成全部轨迹的子轨迹段集合，假设共有lnum条子轨迹段，子轨迹段集合定义为SUB_TRS=｛L₁, L₂, …, L_k, …, L_lnum｝，其中L_k=〈tc_k, tc_k+1〉(1≤k≤lnum)，tc_k, tc_k+1为相邻特征点。

1.2 轨迹上特征点的判定

特征点定义为轨迹中行为变化相对明显的点，能很好地对轨迹结构进行描述，同时应具有一定简缩性及准确性。图像算法中通过检测曲线凹凸分界点来提取曲线边缘信息的方法可以应用到轨迹特征点的提取中。曲线的凹凸分界点指改变曲线向上或向下方向的点。过曲线上任意一点 x 可以作一条关于该条曲线的切线L，切线的方向与曲线上 x 点的方向相同，此时切线L在切点 x 附近的部分“最接近”曲线，由于曲线是其所有切线的包络线，因此在较小范围内，曲线上的点都处于切线L的同一侧。假设给定曲线是由彼此很接近的点集组成，曲线的切线可以由曲线上彼此紧挨着的两个点构成的正向直线代替。据此，推广到轨迹应用中，假设某条轨迹上存在着两个相邻样本点x和x′，且这两个样本点之间的距离非常近，利用无限逼近思想，过这两点可以作一条关于该条轨迹的切线，通过判断样本点与该点相邻轨迹段(即切线)的变化趋势来确定该点是否为特征点。相关定义^[5]如下：

定义4 设平面上两点的坐标分别为P₁(x₁, y₁)和P₂(x₂, y₂)，x₁≠x₂，连接P₁, P₂点得到一条有向直线，称这条有向直线为正向直线，而对应的直线方程:

$\mathit{\boldsymbol{F}}:\left( {{y}_{2}}-{{y}_{1}} \right)(x-{{x}_{1}})+\left( {{x}_{2}}-{{x}_{1}} \right)\left( y-{{y}_{1}} \right)=0$

(1)

为关于P₁, P₂的正向直线方程F。

定义5 假设平面上存在一条正向直线L，根据L将平面上的点分为两类，处于正向直线L顺时针一侧的点称为关于L的内点，而处于正向直线L逆时针一侧的点称为关于L的外点。

定理1 记平面上两点P₁, P₂的正向直线方程F₁₂的左端表达式记为：

${{F}_{12}}\left( x,y \right)=\left( {{y}_{2}}-{{y}_{1}} \right)(x-{{x}_{1}})+\left( {{x}_{2}}-{{x}_{1}} \right)\left( y-{{y}_{1}} \right)$

(2)

对于不在直线L上任一点P₀(x₀, y₀)，则有：

1) 若F₁₂(x₀, y₀)＜0，则点P₀是关于正向直线L的内点；

2) 若F₁₂(x₀, y₀)>0，则点P₀是关于正向直线L的外点。

定理2 设P₁(x₁, y₁)，P₂(x₂, y₂)，P₃(x₃, y₃)，P₄(x₄, y₄)(x₁<x₂<x₃<x₄)，是轨迹上相继的4个点，则：

1) 连接P₁(x₁, y₁)，P₂(x₂, y₂)点得到正向直线L₁的方程为F₁₂(x, y)=0，计算F₁₂(x₃，y₃)的值；

2) 连接P₂(x₂, y₂)，P₃(x₃, y₃)点得到正向直线L₂的方程为F₂₃(x, y)=0，计算F₂₃(x₄, y₄)的值；

3) 若F₁₂(x₃, y₃)*F₂₃(x₄, y₄)＜0，说明在P₃(x₃, y₃)点处轨迹的方向有所改变，认为P₃(x₃, y₃)点是特征点，否则P₃(x₃, y₃)不是特征点。

据此得到轨迹特征点提取方法，即由样本起点P₁P₂构成一条关于轨迹的切线L₁，求P₃点关于这条切线的内点或外点；由P₂P₃点构成另外一条关于轨迹的切线L₂，求P₄点关于这条切线的内点或外点；若P₃和P₄关于各自相邻切线发生变化，即P₃P₄不同为内点或不同为外点，那么认为P₃点是轨迹在P₁P₂P₃P₄处发生变化的点，并添加到特征点集合中。循环判断，直到最后一个样本点。

伪代码如下所示。其中步骤1先将轨迹开始的点P₁和P₂添加到特征点集合中，从第三个点开始判断是否为轨迹特征点，如步骤2到步骤6所示，满足条件的轨迹点将被添加到特征点集合中，不满足条件将依次循环判断下一样本点直到轨迹结束点。步骤7将轨迹结束点添加到特征点集合中，至此完成特征点提取。步骤8至步骤17为轨迹特征点识别方法Calculate_cp()来判断current点是否为特征点，并返回逻辑值。最终利用生成的特征点集合形成该条轨迹的子轨迹段集合。

算法1 Characteristic Points Extraction。

输入：一条轨迹集合TR_i=P₁P₂…P_itnum

输出：该条轨迹特征点集合CHA_TR_i={tc₁, tc₂, …, tc_icnum}

步骤1 将P₁P₂点加入到集合CHA_TR_i中； first=P₁；

步骤2 start=first;current:=start.next.next；

步骤3 while(current.next不为空)，转至步骤8调用方法Calculate_cp计算current是否为特征点；

步骤4 若current是特征点，将current点加入到特征点集合CHA_TR_i中，否则不添加；

步骤5 start=start.next;current=current.next；

步骤6 循环判断结束；

步骤7 将轨迹结束点添加到集合CHA_TR_i中。

//Calculate_cp(start, start.next, current, current.next)

步骤8 初始化cpflag=false；

步骤9 x1:=start.x;y1:=start.y；

步骤10 x2:=start.next.x;y2:=start.next.y；

步骤11 x3:=current.x;y3:=current.y;

步骤12 x4:=current.next.x;y4:=current.next.y；

步骤13 s1:=(x2-x1) *(y3-y1) +(y1-y2) *(x3-x1) ；

步骤14 s2=(x3-x2) *(y4-y2) +(y2-y3) *(x4-x2) ；

步骤15 s3=(x4-x2) *(y3-y2) +(y2-y4) *(x3-x2) ；

步骤16 if(s1*s2) ＜=0) cpflag:=true；

步骤17 return cpflag；

2 轨迹聚类算法 2.1 时空相似性度量

基于时空模式的轨迹聚类顾名思义就是要在时间和空间两种模式下对轨迹数据进行聚类，无论忽略哪一种属性，其聚类结果的准确性及现实意义都将大打折扣，那么针对目前大多数轨迹聚类算法着重考虑轨迹空间特征相似性度量忽略时间特征相似性度量的问题，本文提出一种兼顾时间与空间属性相似性的相似性度量方法，即在空间距离的基础上增加时间相似性距离判断。

时空模式下轨迹的空间相似性度量使用模式识别中常用的空间距离度量方法^{[3, 6]}，包括平行距离、垂直距离和角度距离。假设两条轨迹段分别为L_i(s_i, e_i)以及L_j(s_j, e_j)，其中s_i、e_i、s_j、e_j分别表示L_i和L_j的开始及结束点，并且L_i≥L_j，由短轨迹段L_j两端点s_j、e_j向长轨迹段L_i两端点s_i, e_i垂直投影，得到L_j在轨迹段L_i上的垂直点ps和pe，如图 1所示。

定义6 垂直距离。L_i和L_j之间的垂直距离定义如下：

${{d}_{\bot }}({{\mathit{\boldsymbol{L}}}_{i}},{{\mathit{\boldsymbol{L}}}_{j}})=\frac{l_{\bot 1}^{2}+l_{\bot 2}^{2}}{{{l}_{\bot 1}}+{{l}_{\bot 2}}}$

(3)

定义7 平行距离。L_i和L_j之间的平行距离定义如下：

${d_{//}}({{\bf{L}}_i},{{\bf{L}}_j}) = {\rm{MIN}}({l_{//1}},{l_{//2}})$

(4)

定义8 角度距离。L_i和L_j之间的角度距离定义如下(其中‖L_j‖为L_j的长度)：

${{d}_{\theta }}({{\mathit{\boldsymbol{L}}}_{i}},{{\mathit{\boldsymbol{L}}}_{j}})=\left\{ \begin{matrix} \left\| {{\mathit{\boldsymbol{L}}}_{j}} \right\|\times \sin (\theta ), & 0{}^\circ \le \theta \le 90{}^\circ \\ \left\| {{\mathit{\boldsymbol{L}}}_{j}} \right\|, & 90{}^\circ <\theta \le 180{}^\circ \\ \end{matrix} \right.$

(5)

综上，L_i与L_j之间的空间距离定义如下：

${{d}_{\bot }}({{\mathit{\boldsymbol{L}}}_{i}},{{\mathit{\boldsymbol{L}}}_{j}})+{{d}_{//}}({{\mathit{\boldsymbol{L}}}_{i}},{{\mathit{\boldsymbol{L}}}_{j}})+{{d}_{\theta }}({{\mathit{\boldsymbol{L}}}_{i}},{{\mathit{\boldsymbol{L}}}_{j}})$

(6)

计算轨迹段之间时间距离时，度量单位将是影响聚类分析结果的重要因素。诸多文献都使用标准分数z_score对数据进行预处理，以消除数据间存在着的单位差异^[7-8]。其中文献[8]提出一种不等长时间序列(Short Time Series, STS)滑窗距离聚类算法，该算法使用z_score对时间序列进行标准化预处理，并设计基于滑窗的不等长时间序列距离计算方法，即设计长度等于较短时间段的滑动窗口，从较长时间段开始点向结束点滑动，得到一系列长度等于短时间序列的子序列集合，选择子序列中距离较短时间序列最近的STS距离作为这两条不等长时间序列的时间距离。然而通过实验证明，该方法具有较高的时间复杂度，计算花费非常大。

本文使用一种简单方法对时间数据进行预处理。首先以格林威治标准时间1987-12-31T00:00:00为时间起始点，此时设时间为1，以秒为单位依此类推，即1987-12-31T00:00:01设为2，最终时间范围为168825628~272245926。其次鉴于数值较大计算不便，本文时间数值除以一天秒数(86400) 对时间进行标准化设置，即时间/86400所得值为最终标准化线性时间序列。应用到轨迹时间序列中，假设较长轨迹的时间段为L_i(ist, iet)，ist和iet分别对应L_i的开始时间和结束时间，较短轨迹的时间段为L_j(jst, jet)，其中jst和jet分别对应L_j的开始时间和结束时间，‖L_i‖≥‖L_j‖。比较两条轨迹的时间段，可得三种情况：L_i完全包含L_j，L_i部分包含L_j以及L_i完全不包含L_j，三种情况如图 2(a)所示，其中‖L_j‖、‖L_j′‖、‖L_j″‖分别表示其长度。

应用到实际情况中，可能存在两条轨迹段时间跨度相差较大的情况，在这里取两条轨迹段之间的最近距离作为两条轨迹段的时间距离。将较短的子轨迹时间段L_j向较长的子轨迹时间段L_i投影，以左投影为例，如图 2(b)所示，将较短时间段结束点jet投影至较长时间段上距离其开始点ist长度为|L_j|的点，即iet′点，其中|iet′-ist|=|L_j|。据此将L_i与L_j之间的时间距离定义如下：

$DT=\left\{ \begin{array}{*{35}{l}} 0, & {{\mathit{\boldsymbol{L}}}_{j}}\subset {{\mathit{\boldsymbol{L}}}_{j}} \\ \sqrt{{{\left| ist-jst \right|}^{2}}+{{\left| ie{t}'-jet \right|}^{2}}}, & 其它 \\ \end{array} \right.$

(7)

考虑到空间距离和时间距离存在着不同数量级的问题，在此使用z-score标准分数对空间及时间距离进行标准化。对变量f的度量值进行如下标准化^[9]。

1) 计算绝对偏差的平均值：

${{S}_{f}}=\frac{1}{n}\left( |{{x}_{1f}}-{{m}_{f}}|+|{{x}_{2f}}-{{m}_{f}}|+\cdots +|{{x}_{nf}}-{{m}_{f}}| \right)$

(8)

其中：x_1f, x_2f, …, x_nf为f的n个度量值，m_f为f的均值:

${{m}_{f}}=\frac{1}{n}\left( {{x}_{1f}}+{{x}_{2f}}+\cdots +{{x}_{nf}} \right)$

(9)

2) 计算标准度量值：

${{Z}_{if}}={{x}_{if}}-{{m}_{f}}/{{s}_{f}}$

(10)

Z_if为标准化后度量值，表示原始数据偏离平均值的程度, 且服从正态分布。设标准化后空间距离和时间距离分别为DS′、DT′，通过设置距离权重w(0≤w≤1) ，可以用来调整对空间距离和时间距离的敏感度，那么轨迹段间时空距离定义为：

$DST=w\times DS'+(1-w)\times DT'$

(11)

2.2 时空轨迹聚类算法

通过观察轨迹分段得到的子轨迹段集，可以发现子轨迹段形状具有不规则性，且集合中含有大量的噪声。由于DBSCAN算法通过对区域密度的连通性分析来聚类，不仅可以发现任意形状的聚类簇，而且在聚类过程中能够最大限度地避免噪声的干扰，因此在此采用该方法对子轨迹段进行聚类。该算法需要设定两个全局参数，即邻域半径eps和minlns。DBSCAN通过筛选子轨迹数据集中每个对象的eps邻域来搜索簇，如果对象p的eps邻域中包含的对象个数大于等于minlns个，则创建一个以p为核心对象的簇。然后，算法迭代地聚类从这些核心对象直接密度可达的所有对象，这个过程中可能涉及一些密度可达簇的合并，当没有新的对象添加到任何簇时，过程结束。下面给出相关定义^[10]。

定义9 eps邻域。以给定对象p为中心，半径为eps的区域称为该对象的eps邻域。

定义10 核心对象。在给定对象p的eps邻域内，如果样本点数大于等于minlns，则称该对象为核心对象。

定义11 直接密度可达。若给定对象q为核心对象，并且对象p在核心对象q的eps邻域内，那么称对象p从对象q出发是直接密度可达。

定义12 密度可达。样本集合D中，如果存在一个对象链p₁，p₂，…，p_n，设p₁=q，p_n=p，若对象链中任意p_i+1是从p_i出发关于eps和minlns直接密度可达的，那么对象p是从对象q出发关于eps和minlns密度可达。

定义13 密度相连。如果样本集合D中存在对象o，使得对象p和q都是从o出发关于eps和minlns密度可达，那么称对象p和o是从对象q出发关于eps和minlns密度相连。

应用到轨迹段聚类中，需要扫描整个轨迹段数据集，找到任意一个核心轨迹段，对该核心轨迹段进行扩充。扩充的方法是寻找从该核心轨迹段出发的所有密度相连的轨迹段。遍历该核心轨迹段的eps邻域内的所有核心轨迹段，寻找与这些轨迹段密度相连的对象，直到没有可以扩充的轨迹段为止。之后重新扫描数据集(不包括之前寻找到的簇中的任何轨迹段)，寻找没有被聚类的核心轨迹段，再重复上面的步骤，对该核心轨迹段进行扩充直到数据集中没有新的核心轨迹段为止。数据集中没有包含在任何簇中的轨迹段可以认为是噪声点。算法时间复杂度的限制主要体现在邻域的计算，在不使用索引的情况下邻域计算需要遍历所有轨迹段，时间复杂度为O(n²), 导致计算效率低；使用索引能够降低时间复杂度，如文献[3-4]都使用R-tree结构进行索引来提高整体计算速度，而Gudmundsson^[11]近期提出利用并行算法来计算子轨迹集群的方法，以改进现有串行方法为基础，利用结合了GPU计算能力的并行技术来实现算法速度的大幅提高。由于其实现代价比较昂贵，在这里提出另外一种方法来降低时间复杂度。参考文献[12]设计交通路网的思想，本文首先根据数据样本集定义m*n的网格空间，并将每一格空间进行序号标记，统计每条轨迹穿过的所有网格。计算某条子轨迹段邻域的时候先以该条轨迹段所在网格为中心，向四周扩展，形成3×3的九宫格搜索区域，最终只需判断所有穿过该九宫格区域的轨迹，示意图如图 3所示。

由此，本文提出一种改进的轨迹聚类算法，伪代码如下：

算法2 Trajectory Segments Clustering(TSC)。

输入：轨迹段集合TS={L₁, L₂, …, L_lnum}, 邻域半径eps, 最小轨迹段数minlns。

输出：聚类结果集合C={C₁, C₂, …, C_clusnum}。

步骤1 初始化clusterId=0, 并将轨迹集合中所有轨迹段初始化为 unclassified。

步骤2 循环判断集合TS中所有轨迹段，若轨迹段L_i的标识为 unclassified 则执行步骤3。

步骤3 根据eps计算L_i的邻域N_ε(L_i)。

步骤4 若|N_ε(L_i)|≥minlns，则为L_i所有邻域轨迹段N_ε(L_i)分配为当前clusterId，并将N_ε(L_i)-{L_i}添加到待扩展序列Q中，执行步骤8对序列Q进行可达近邻扩展；否则将此L_i设置为noise。

步骤5 clusterId++；将此L_i设置为classified。

步骤6 循环判断结束。

步骤7 根据clusterId，输出聚类结果集合C={C₁, C₂, …, C_clusnum}。

步骤8 while(Q不为空)，选择Q中任意Q_i。

步骤9 根据eps计算Q_i的邻域N_ε(Q_i)。

步骤10 若N_ε(Q_i)≥minlns, 执行步骤11；否则将Q_i从序列Q中移除。

步骤11 对每一个X(其中X∈N_ε(Q_i))，若X为noise, 将当前clusterId分配给X；否则将X添加至序列Q中。

步骤12 循环判断，直到Q为空。

3 实验结果

为测试所提算法的有效性，使用两个真实数据集进行测试：一是Starky收集的动物迁徙数据集(http://www.fs.fed.us/pnw/starkey/data/tables/)，提取该数据集1993年麋鹿的运动轨迹(简称Elk93) ，并从传感定位数据源中提取经纬度坐标以及定位时间信息参与实验分析；Elk93由47204个采样点构成33条轨迹。二是大西洋飓风数据(http://weather.unisys.com/hurricane/atlantic/)，选取经度纬度、采样时间属性参与实验；实验中选取1950—2009年之间的数据，经整理数据集由19750个采样点组成639条轨迹。

为直观体现在轨迹提取特征点阶段所提方法的合理性，实验中选取1950年大西洋飓风数据中的三条轨迹进行测试(共124个样本点)，如图 4所示，图中空心圆为样本点，整条轨迹由彼此紧挨着的相邻样本点构成。

在实验中，使用三种提取特征点的方法进行实验比较，分别为文献[3]提出的一种近似MDL方法和文献[4]提出的利用相邻样本点转角角度变化来提取特征点的方法(其中转角阈值设置为45°)，以及本文提出的利用曲线边缘检测提取特征点的方法，提取结果如图 5所示。以拥有32个样本点的2号轨迹为例，MDL算法提取了11个特征点(+标注)，简缩率约65%，而“转角”方法提取3个特征点(□标注)，简缩率90.6%，曲线边缘检测法提取特征点4个，简缩率约87.5%。

三种特征点提取方法比较结果如表 1所示。三条轨迹共有124个样本点，使用MDL方法提取特征点数为30个，在精确性方面表现良好，但在简缩性方面有些不足；由于轨迹相邻样本点转角变化不大导致根据转角角度变化来提取特征点的方法表现欠佳，提取特征点数为7个，特别是1号和3号轨迹只提取到了轨迹的开始点和结束点，不能对轨迹进行准确描述。本文利用曲线边缘检测方法的原理，提取出的特征点数为14个，在保持原轨迹精确性和子轨迹简洁性之间找到最优的平衡，表明了所提方法的合理性。

为了比较不同方法轨迹聚类的质量，定义如下聚类质量评价公式^[3]：

$QMeasure=\sum\limits_{i=1}^{nu{{m}_{clus}}}{\left( \frac{1}{2|{{C}_{i}}|}\sum\limits_{x{{\in }_{Ci}}}{\sum\limits_{y{{\in }_{Ci}}}{dist{{(x,y)}^{2}}}} \right)}$

(12)

在时空模式下，通过计算所有聚类簇中相互轨迹段间的平方误差之和来衡量轨迹聚类质量。QMeasure越小表示聚类簇内轨迹段密度越大，聚类结果质量越高。其中dist(x, y)为聚类簇中两条轨迹段间的距离，如式(11) 所示。为取得较好结果，在不同权重下分别对两种数据进行实验。图 6为eps=9、minlns=30时，对飓风数据在不同权重下进行聚类所得聚类结果质量分析图。从图中可以看到，w=0.7时聚类质量达到最优，整体对权重设置具有一定的敏感性。

图 7为eps=8、minlns=28时，对ELK93进行不同权重下聚类质量的分析比较。聚类质量在w=0.5达到最优，然而由于ELK93数据采样自1993年春季到秋季期间麋鹿的运动轨迹，时间跨度较小，故而对权重设置敏感性较低。

在此，实现了轨迹聚类(TRAjectory CLUStering, TRACLUS)算法和根据转角进行轨迹划分并由DBSCAN完成聚类的算法(简称“转角”算法)，以及本文中提出的改进算法。图 8为eps=9时对大西洋飓风数据的聚类结果质量分析图，此时w为0.7。

从图 8中可以看出，使用TRACLUS算法对飓风数据进行聚类，在minlns=31时达到质量最优，但由于TRACLUS在进行相似性度量时只考虑了空间特征，导致整体性能低于本文中提出算法。“转角”算法中由于根据相邻样本点变化角度来划分轨迹段结果质量低下导致QMeasure值较前两种方法偏高，为解决这种情况，可事先设定转角累加阈值，相邻样本点转角低于转角阈值时将之累加直至达到累加阈值，这时将样本点添加为特征点同时累加阈值清零，但由于需要增加参数值设置，且对结果影响较大，因此本文没有采用这种方法。使用本文算法聚类质量在minlns=30时达到最优，参数设置对聚类质量有一定的影响，但聚类质量值整体较前两种算法偏低。

图 9显示了eps=8时对ELK93数据进行聚类的聚类质量图，此时w为0.5。从图中可以看出使用TRACLUS算法对ELK93数据进行聚类，聚类质量在minlns=28时达到最优，而minlns过大或过小都会对聚类质量造成较大影响。“转角”算法结果与TRACLUS类似，在minlns=28时达到最优，总体呈“倒U型”，但聚类质量低于TRACLUS。使用本文算法进行聚类，聚类质量值从minlns=25开始持续下降，在minlns=29时达到最优，因为同时考虑了空间及时间相似性，所以聚类质量总体较高。

以ELK93 为例，同样使用文献[3]中TRACLUS算法以及本文提出改进聚类算法进行比较，比对结果如图 10所示，其中(a)为使用TRACLUS算法进行聚类得到聚类簇代表轨迹，(b)为应用本文改进算法得到的聚类结果，从图中可以看出使用TRACLUS算法在只考虑空间相似性的情况下，生成聚类簇相比本文算法较大，簇内密度小。观察图 10(a)最下方的一条聚类簇代表轨迹，相应本文算法将之生成为三个相对较小聚类簇，说明本文在同时度量轨迹数据时空相似性的情况下，生成聚类簇更加精细，更有使用价值。

4 结语

对轨迹进行划分再聚类的方法能够很好地发现存在于复杂时空轨迹中的局部相似性，然而其聚类结果对子轨迹段划分质量具有很强的依赖性，而轨迹特征点是影响子轨迹段划分质量的关键因素，同时相似性度量也对聚类结果具有重要影响。针对相似性度量中忽略时间特征的问题，提出了一种基于时空特征的轨迹数据聚类方法。该方法借鉴图像算法中曲线边缘信息的检测方法，通过计算轨迹上某一点与该点相邻轨迹段的变化趋势来判定该点是否为轨迹的特征点。实验表明使用该方法所提取的特征点具有较高的简洁性，并且根据特征点生成的子轨迹段集能够准确地描述轨迹结构，降低了轨迹划分对聚类质量的影响。在对子轨迹段集进行聚类的过程中，兼顾轨迹的空间与时间特征，并通过设置权重调整时空特征敏感度，实验表明基于时空模式的轨迹数据聚类方法能够发现移动对象内在结构及其动态规律。接下来的工作是进一步提高聚类效率。