轨迹数据在诸多实际领域都有十分重要的作用，如智能交通等[1]。Sermanet等[2]使用聚类算法对行人轨迹数据进行分析。Le[3]使用多尺度模型对轨迹数据进行无监督聚类。为了获取数据的模型或发现数据中隐藏的深层信息，聚类是分析模型常用且有效的途径[4-5]。目前聚类算法已经有了大量的研究，常见的聚类算法有k-means[6]、BIRCH (Balanced Iterative Reducing and Clustering using Hierarchies)[7]、具有噪声的基于密度的聚类方法 (Density-Based Spatial Clustering of Applications with Noise, DBSCAN)[8]、OPTICS (Ordering Points To Identify the Clustering Structure)[9]以及STING (STatistical INformation Grid)[10]等。在上述聚类算法中，k-means算法存在一个棘手的问题，即如何在聚类前确定数据中类的个数，而且异常数据对聚类效果的影响较大。BIRCH算法对数据分布的形状有一定的要求。DBSCAN算法的邻域半径参数需要人工给出。OPTICS算法虽然对DBSCAN算法进行了改进，但其内部复杂处理过程使得实际运行速度远远低于DBSCAN。STING算法形成的类簇的外边界不确定性较大。这些聚类算法难以直接应用在复杂轨迹数据上，主要因为轨迹数据有以下三点特殊性：每条轨迹都是一个单独的元素，其内部拥有不同的采样点个数；在同一场景的不同区域里，轨迹间的相似性不尽相同；正常轨迹与异常轨迹往往会存在交叉重叠等可能性。
Vlachos等[11]使用了小波变换改变轨迹的表现形式来提升聚类算法效果。在图像领域，基于块匹配的三维滤波算法 (Block-Matching and 3D filtering, BM3D) 是目前较为领先的图像滤波算法[12]，该算法首次利用了相似结构的图像块的一致性使得滤波效果得以提升。受以上核心思想的启发，本文提供了一种对轨迹做摘要处理的算法，通过对轨迹重采样以及非局部去噪，寻找更有效的轨迹数据聚类和异常检测算法。
1 轨迹摘要算法框架本文算法框架主要是由两部分构成，即重采样和非局部去噪。下面的章节将详细介绍这两个过程。
1.1 重采样通常不同的轨迹有不同个数的采样点，难免包含一些噪声采样点，会对后续距离计算造成不准确的影响，且后续的硬阈值滤波过程需要使用的是等长 (即等采样点个数，下文中所有等长均为此含义) 轨迹。本文使用了重采样技术，使得轨迹得以平滑处理，将噪声采样点的影响降低，并使每一条轨迹等长。每一条轨迹的采样点都会被重采样为相同的间隔。为了使重采样过程对原始轨迹数据信息量的破坏降到最小，需要针对一个轨迹数据集找到一个最优采样点数moptimal，尽可能保留原始轨迹的数据特征。
在本文中，重采样所处理的轨迹数据是二维坐标点序列，不妨令Traj={p1, p2, …, pm}表示一条具有m个采样点轨迹，此处的pi=(xi, yi) 是二维坐标。本文将二维坐标序列中的时间序号i作为自变量，将x坐标和y坐标分别作为因变量。此时，对 (i, xi) 序列做多项式拟合，使用最小二乘法来优化多项式的参数。同样的，也对 (i, yi) 序列拟合，根据实验结果表明，通常多项式拟合中的最高项次数设置为8比较合适。分别将x曲线和y曲线划分为moptimal-1段，可以计算出moptimal个x坐标以及moptimal个y坐标，由此可以得到一条具有moptimal个采样点的重采样轨迹。
考虑到重采样后需要保留原始轨迹的整体形状信息，为获取最优值moptimal，使用重采样前后的JSD (Jensen-Shannon Divergence) 距离作为衡量标准，目的是最小化JSD距离。因此在这里将轨迹采样点的夹角分布作为轨迹形状的重要特征，以考量轨迹重采样前后的信息变化量。
轨迹的夹角信息{θi, i=1, 2, …, m-2}如图 1所示。
|
图 1 轨迹点和角度说明 Figure 1 View of trajectory points and angles |
使用具有B个bins的角度直方图建立概率分布P(x)，夹角的角度值范围从min{θi}至max{θi}，同时设置B=12。通过实验，发现moptimal的最优值对B的变化并不敏感。对重采样后的轨迹Traj′，可以使用同样的方法取得Q(x)，Traj和Traj′的JSD距离可以通过式 (1) 计算:
| $JSD(P||Q)=\frac{1}{2}D(P||M)+\frac{1}{2}D(Q||M)$ | (1) |
其中:
|
图 2 JSD曲线 Figure 2 JSD curve |
在轨迹重采样之后，本文使用非局部去噪方法，通过多次迭代获得多尺度摘要。该方法包含冗余匹配、硬阈值处理以及合并3个组成部分。异常轨迹会在某一次迭代中被剔除，且不会参与后续的迭代过程。
为了便于理解，定义TRk作为第k次迭代的输出轨迹集合，同时也作为第k+1次迭代的输入轨迹集合。TR0是原始轨迹重采样后的轨迹集合。迭代的终止条件是两次相邻迭代的输出轨迹结果几乎不改变，即TRk-1≈TRk。
1.2.1 冗余匹配在这一步骤中，相似的轨迹将被互相匹配并建立轨迹组，与此同时，异常轨迹也可能在相似性匹配的过程中被检测和剔除。每一条轨迹会作为一个参考轨迹，匹配与其相似性较高的其他轨迹。假定当前参考轨迹为Traj={p1, p2, …, pm}，则其相似轨迹组G可以通过式 (2) 计算得到:
| $\mathit{\boldsymbol{G}}=\left\{ Tra{{j}_{x}}|dis\mathit{tan}ce\left( Tra{{j}_{x}},Traj \right)\le \xi \right\}$ | (2) |
由于轨迹是等长的，在式 (2) 中，distance函数为轨迹间欧氏距离，即将两条轨迹的所有对应点的欧氏距离求均值。在匹配过程中，冗余的体现之处在于一条轨迹的信息可能会出现在多个轨迹组当中，且在各个组内的处理是独立的。ξ是距离的阈值，可以自适应地选取，在后续章节本文将详细阐述自适应的方法。异常轨迹所在的分组可以较为容易地和其他轨迹组区分，通常如果轨迹组内的轨迹数量不足3条，则认为该条参考轨迹是异常的。
1.2.2 硬阈值处理在本过程中，需对所有轨迹组进行小波硬阈值[13]处理。假定Traj为一条参考轨迹，拥有g个相似的邻居轨迹，则轨迹组G信号形式为:
| $\mathit{\boldsymbol{G}}=\left[ \begin{matrix} Tra{{j}_{1}} \\ Tra{{j}_{2}} \\ \vdots \\ Tra{{j}_{g}} \\ \end{matrix} \right]=\left[ \begin{matrix} {{p}_{11}} & {{p}_{12}} & \ldots & {{p}_{1m}} \\ {{p}_{21}} & {{p}_{22}} & \ldots & {{p}_{2m}} \\ \vdots & \vdots & {} & \vdots \\ {{p}_{g1}} & {{p}_{g2}} & \ldots & {{p}_{gm}} \\ \end{matrix} \right]={{\left[ \begin{matrix} {{s}_{1}} \\ {{s}_{2}} \\ \vdots \\ {{s}_{m}} \\ \end{matrix} \right]}^{\rm{T}}}$ | (3) |
此处si为轨迹组的第i个信号，由组内所有轨迹的第i个采样点组成，之后对si作小波硬阈值处理。首先，对si做小波变换，将其转换为小波系数；其次，对小波系数作硬阈值处理，使得数据更为平滑；最后，去噪之后的信号会通过小波反变换还原到原始空间，整个过程如式 (4) :
| ${{s}_{i}}=wavele{{t}^{-1}}\left( threshold\left( \text{ }wavelet\left( {{s}_{i}} \right) \right) \right)$ | (4) |
其中:wavelet(·) 表示小波变换，threshold(·) 表示硬阈值滤波，wavelet-1表示小波反变换。
1.2.3 合并在之前的过程中，对每条轨迹都作了冗余匹配，因而在硬阈值处理之后，不同的轨迹组中的轨迹可能来源于本次迭代分组前的同一条轨迹Trajx。将这些轨迹合并，作为轨迹Trajx本次迭代的结果。合并过程即为将这些轨迹相对应采样点求加权均值，权重为轨迹所在分组的轨迹数量。使用本次迭代的轨迹结果作为下一次迭代的输入数据，直到满足迭代的终止条件。
1.3 参数的自适应选取对于某个轨迹数据集，距离阈值ξ可以通过数据集的所有轨迹自适应地获取。首先，对于每一条轨迹，可以获得该条轨迹到其他所有轨迹的欧氏距离，并将这些距离按升序排序。其次，将所有轨迹的第x小距离相加，并求平均值，结果如图 3(a)所示。该图像是单调递增的，且增长速度在不同的位置也不同。最后对该函数求二阶导，即求出增长速度变化量的曲线，如图 3(b)所示，图中的最高点表明距离变化最大之处，此处x=153。将阈值ξ设置为图 3(a)中x=153时y坐标的值。这是因为，小于ξ的部分代表类簇内部的距离，大于ξ的部分代表类簇之间的距离，在此处会出现距离增大速度十分陡峭的现象。阈值ξ用作轨迹相似性的判定是合适的。
|
图 3 阈值ξ的选取 Figure 3 Selection of threshold ξ |
为了更清晰地体现本文提出的轨迹摘要算法的特点，本文对两个数据集进行了可视化，以便分析多次迭代的轨迹结果。
2.1 交叉数据集的分析图 4(a)是交叉路口数据集[1]的原始轨迹，数据集中有轨迹262条，轨迹的采样点个数区间为[8, 22]。迭代结果分别为图 4(b)~图 4(f)。从宏观上看，轨迹数据从杂乱重叠的原始数据，被压缩为了精简的7类轨迹，也可以成为7类模式。迭代在第5次之后停止，因为这次迭代的结果和第4次的迭代结果并无差异，满足了迭代的终止条件。从图中可以看出，本身较为复杂的轨迹数据通过多次迭代形成了多尺度的结果，例如图 4(a)~图 4(b)的过程，轨迹被大幅度地规整，并且其中有12条异常轨迹被检测出，因而在第1次迭代的结果中被剔除了；在图 4(b)~图 4(c)的第2次迭代中，有一条从左向右贯穿的轨迹，由于其右侧明显的拐角变化 (箭头指向位置)，因而被检测为异常轨迹。图 4(g)是DBSCAN的聚类结果，本文使用不同的灰度值区分不同的类，位置明显不同的轨迹也属于不同的类。可以看到，DBSCAN将整个轨迹数据集聚为5类。如图 4(h)，从DBSCAN的聚类结果的每个类中选出一条能代表该类的中心轨迹，作为数据集的摘要。在图 4(f)中，原始轨迹在1处和2处的轨迹数量并不多，但也有一定的规模，不能将其简单地视为异常，本文算法最终保留下了这些信息，而在DBSCAN的结果中这些轨迹则被视为了异常。
|
图 4 交叉路口数据集摘要结果 Figure 4 Abstraction of cross dataset |
图 5(a)是实验室数据集[1]的原始轨迹，轨迹条数为122，采样点个数的范围是[39, 624]。各次迭代的轨迹摘要结果分别为图 5(b)~图 5(f)。在图 5(f)中，可以清晰地看到图中7条轨迹，这7条中的每一条轨迹实际上是数十条相同的轨迹重合的结果。图 5(g)和图 5(f)分别是DBSCAN的聚类结果和聚类中心。
|
图 5 实验室数据集摘要结果 Figure 5 Abstraction of lab dataset |
本文提出的轨迹摘要算法，不仅能提供简洁的可视化结果，也能通过多次迭代对轨迹数据进行异常检测。为评估本文算法的异常检测效果，本文引入了召回率 (Recall) 和准确率 (Precision) 两种指标，与DBSCAN算法作对比。之所以引入这两个指标，是因为在异常检测算法的结果中，可能有一些异常的数据没有被检测到，也可能把一些正常的数据视为异常。Recall指标的意义是衡量异常的数据是否被检测到的情况，Precision指标的意义是衡量把正常的数据是否视为异常的情况。假定A表示算法找到的异常轨迹中确实是真实的异常轨迹的数量，B表示数据中的真实异常轨迹的总数量，C表示算法所认为的异常轨迹的数量。Recall可以表示为A/B，即算法找出的真正的异常轨迹条数与数据集中的真实异常轨迹条数的比值；Precision可以表示为A/C，即算法找出的真正的异常轨迹条数与算法所认为的异常轨迹数量的比值。为客观评价，两种算法均采用自适应方法选取距离阈值。本文所有实验在Intel Core i5-3570 CPU 3.4 GHz，内存4 GB上运行。具体实验结果见表 1，这两种指标的值越大表明异常检测的能力越好。
|
表 1 2种算法的性能指标对比 Table 1 Performance comparison of two algorithms |
在表 1中，本文算法并不是在每一项指标中都完全优于DBSCAN算法，因为Recall和Precision两种指标的侧重点不同，且它们之间存在制约关系，所以可能出现此消彼长的情况。为了对比算法异常检测的综合性能，可以使F1指标 (F1-measure), 该指标由式 (5) 计算:
| $F1=2PR/\left( P+R \right)$ | (5) |
其中：R表示Recall，P表示Precision，该指标的对比结果见表 1。从表 1可以看出，本文算法在各个数据集上异常检测的综合性能基本优于DBSCAN算法。
4 结语为了对复杂轨迹数据聚类和异常检测，本文提出了一种轨迹摘要算法。通过重采样，将轨迹等长化，为后续算法提供了可用的数据，同时对可视化效果也有所提升; 采用自适应阈值，使轨迹相似性的计算更为准确; 利用非局部去噪，将整个轨迹数据集拆分为各个小部分，以便独立处理，提升了异常检测效果，算法运行效率也得以提升。在后续研究中，主要工作是考虑将较长轨迹分段处理，同时改进当前异常检测的方法，以及尝试提升轨迹摘要可视化效果的方法。
| [1] | MORRIS B T, TRIVEDI M M. Understanding vehicular traffic behavior from video:a survey of unsupervised approaches[J]. Journal of Electronic Imaging, 2013, 22(4): 041113. doi: 10.1117/1.JEI.22.4.041113 |
| [2] | SERMANET P, KAVUKCUOGLU K, CHINTALA S, et al. Pedestrian detection with unsupervised multi-stage feature learning[C]//Proceedings of the 2013 IEEE Conference on Computer Vision and Pattern Recognition. Piscataway, NJ:IEEE, 2013:3626-3633. |
| [3] | LE Q V. Building high-level features using large scale unsupervised learning[C]//Proceedings of the 2013 IEEE International Conference on Acoustics, Speech and Signal Processing. Piscataway, NJ:IEEE, 2013:8595-8598. |
| [4] | LAXHAMMAR R, FALKMAN G. Online learning and sequential anomaly detection in trajectories[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2014, 36(6): 1158-1173. doi: 10.1109/TPAMI.2013.172 |
| [5] | MORRIS B T, TRIVEDI M M. A survey of vision-based trajectory learning and analysis for surveillance[J]. IEEE Transactions on Circuits and Systems for Video Technology, 2008, 18(8): 1114-1127. doi: 10.1109/TCSVT.2008.927109 |
| [6] | KHAN S S, AHMAD A. Cluster center initialization algorithm for K-means clustering[J]. Pattern Recognition Letters, 2004, 25(11): 1293-1302. doi: 10.1016/j.patrec.2004.04.007 |
| [7] | ZHANG T, RAMAKRISHNAN R, LIVNY M. BIRCH:an efficient data clustering method for very large databases[C]//Proceedings of the 1996 ACM SIGMOD International Conference on Management of Data. New York:ACM, 1996:83-94. |
| [8] | ESTER M, KRIEGEL H P, SANDER J, et al. A density-based algorithm for discovering clusters in large spatial databases with noise[J]. Knowledge Discovery and Data Mining, 1996, 96(34): 226-231. |
| [9] | ANKERST M, BREUNIG M, KRIEGEL H P. Ordering points to identify the clustering structure[EB/OL].[2016-06-20].http://www.nwadmin.de/folien-optics.pdf. |
| [10] | WANG W, YANG J, MUNTZ R. STING:a statistical information grid approach to spatial data mining[C]//Proceedings of the 23rd International Conference on Very Large Data Bases. San Francisco, CA:Morgan Kaufmann Publishers, 1997, 97:186-195. |
| [11] | VLACHOS M, LIN J, KEOGH E, et al. A wavelet-based anytime algorithm for k-means clustering of time series[EB/OL].[2016-06-20].http://cs.gmu.edu/~jessica/publications/ikmeans_sdm_workshop03.pdf. |
| [12] | DABOV K, FOI A, KATKOVNIK V, et al. Image denoising by sparse 3-D transform-domain collaborative filtering[J]. IEEE Transactions on Image Processing, 2007, 16(8): 2080-2095. doi: 10.1109/TIP.2007.901238 |
| [13] | JOHNSTONE I M, SILVERMAN B W. Wavelet threshold estimators for data with correlated noise[J]. Journal of the Royal Statistical Society:Series B, 1997, 59(2): 319-351. doi: 10.1111/rssb.1997.59.issue-2 |
| [14] | ANJUM N, CAVALLARO A. Multifeature object trajectory clustering for video analysis[J]. IEEE Transactions on Circuits and Systems for Video Technology, 2008, 18(11): 1555-1564. doi: 10.1109/TCSVT.2008.2005603 |