时空轨迹多维特征融合的行为规律挖掘算法

2023-04-19姜乔文刘瑜谭大宁孙顺董凯

航空学报 2023年5期

姜乔文，刘瑜，*，谭大宁，孙顺，董凯，2

1．海军航空大学信息融合研究所，烟台 264000

2．中国电子科学研究院，北京 100041

随着雷达、卫星和无线传感器网络等目标定位技术的迅猛发展，已经基本形成了覆盖天基、空基、岸基、海上和水下等多层次的传感器预警监测体系［1-2］，各种被监视目标的行为得到了准确的检测、跟踪和关联，最终以时空轨迹的形式记录下来。海量的时空轨迹数据在预警监视系统中不断地生成和积累，其中蕴含着大量规律性的信息和知识。当前目标监视数据库已经初具规模，但是由于各型传感器的探测感知能力有限，在实际监视过程中获得的时空轨迹往往存在采样率不稳定、数据多源异质等问题，这对于多源探测数据的融合处理带来了巨大的挑战［3］。因此，基于时空轨迹数据多维特征挖掘目标行为规律，分析和预测目标行为意图，对于提升信息质量、辅助指挥决策、实现智能化态势感知具有重要意义。

在时空轨迹数据挖掘领域，通过对目标海量的历史轨迹数据进行相似性度量和聚类分析，可以发现目标在运动过程中呈现出的相似性行为模式。近年来，国内外大量专家学者在轨迹相似性度量［4-8］和轨迹聚类［9-14］方面进行了大量有价值的研究。周星星等［15］对经典的时空轨迹相似性度量方法进行了总结和分类，分析了各自优缺点和应用场景。Liao［16］总结了时间序列聚类的概念、距离度量、聚类算法、评价方法以及应用领域。Zheng［17］系统地综述了现有的轨迹数据挖掘技术，构建了从轨迹数据预处理到轨迹相似性度量再到轨迹聚类的目标行为模式挖掘框架，并提供了一些公开的轨迹数据集。虽然目前时空轨迹行为规律挖掘方法的研究已经取得了很大进展，但是现有研究大部分仅考虑轨迹的空间位置信息，忽略了时间、航向和速度等多维信息，因而难以区分空间位置相似但运动速度和方向不同的行为。Pan 等［18］利用轨迹的位置、航向和速度信息提出了基于多维轨迹特征的行为规律挖掘算法，但是该算法通过权重系数融合量纲不同的特征信息，存在参数设置困难、稳定性不强等缺点；除此之外，该算法容易受到轨迹行为分布密度的影响，易造成虚警和漏警。

针对现有轨迹行为规律挖掘方法的局限性，本文提出了一种时空轨迹多维特征融合的行为规律挖掘算法，通过轨迹的时间和空间信息进行多维特征融合相似性度量，基于k最近邻距离和决策图进行轨迹聚类，实现未知场景下对目标行为规律的无监督挖掘。

1 时空轨迹

时空轨迹数据是由若干个目标数据点按时间顺序组成的序列。假设预警监视系统中积累的所有目标时空轨迹数据为

式中：TR 表示时空轨迹数据集；i为航迹编号；n为航迹总条数；tri表示数据集中的第i条时空轨迹，即

其中：pi，k表示第i条时空轨迹中的第k个目标信息数据点；m表示时空轨迹tri中信息点总数。由于各种传感器的采样率存在差异，加上偶尔的探测数据缺失，因此在不同的时空轨迹中，即使是运动模式相似的轨迹，信息点总数通常也并不相同。

pi，k是一个含有多维信息的目标数据点，其中至少包含探测时刻ti，k和目标位置li，k这2个重要信息，ti，k＜ti，k+1，li，k为二维平面上或三维空间中的一个点。除此之外，根据传感器的探测感知能力差异，其中还可能包含在该时刻目标的航向θi，k、速度vi，k、加速度ai，k等特征，可以将其表达为一个多维向量［18］：

对于仅含有ti，k和li，k信息的时空轨迹，既可以把它理解为带有时间标记的空间序列，也可以理解为带有空间标记的时间序列［19］。Hung 等［20］将时空轨迹视为一种时空棱柱结构，例如一个在二维平面上运动的目标，X轴和Y轴构成的平面为其位置空间，li，k中的位置信息为(xi，k，yi，k)，T轴表示其时间信息，那么它的时空轨迹tri在XY平面上的投影即为位置信息点li，k构成的轨迹曲线，如图1 所示。

图1 时空轨迹示意图Fig.1 Schematic diagram of space-time trajectory

2 时空Hausdorff 距离

2.1 时空轨迹相似性度量

在时空轨迹数据挖掘中，相似性度量是一个非常关键的过程，决定着后续轨迹聚类结果的好坏。轨迹间的相似性度量方法是通过衡量不同轨迹在空间上的远近程度来判断是否相似，由于目标行为规律挖掘的应用场景和目的需求的不同，相似性度量算法对于轨迹间的匹配度有着不同的解释。目前应用广泛的方法主要有基于欧式距离（Euclidean Distance， ED）的相似性度量［4］、编辑距离［5］（Edit Distance on Real sequence， EDR）、豪斯多夫距离［6］（Hausdorff Distance，HD）、最长公共子序列［7］（Longest Common Sub-Sequence， LCSS）、动态时间归整［8］（Dynamic Time Warping， DTW）等以及其改进算法。Agrawal 和Swami［21］运用欧氏距离来计算轨迹对应点的相似程度，其优点在于计算简单，复杂度只有O（n），但要求轨迹的采样率和点数量一致，在实际应用中几乎无法满足条件。魏龙翔等［22］结合Hausdorff 距离和LCSS 同时利用轨迹的位置信息和方向信息进行相似性度量。Li等［23］针对DTW 算法在距离度量使存在过度拉伸和压缩的问题，通过引入惩罚函数提出了一种自适应约束DTW 算法，可以自适应的调整时空轨迹间的对应关系。Lee 等［24］改进了欧氏距离的相似性度量方法，提出了水平、垂直以及角距离等多个维度融合的Hausdorff 距离度量算法。

不同于匀速直线运动、加速直线运动、转向等简单的、可以模型化的目标局部行为模式，本文的行为规律是指目标根据其本身特性和运动目的在某一时段、某一区域呈现出的全局运动特点。例如在空管交通中，军用飞行器和民用飞机由于性质和运动目的不同，在同一区域的全局运动轨迹和运动方式呈现出不同的模式，这就需要综合利用时空轨迹的多维特征来进行相似性度量。然而目前大部分度量算法仅考虑轨迹的空间位置信息，忽略了时间、航向和速度等多维信息，导致后续的聚类算法无法挖掘出空间位置相似但运动速度和方向不同的行为，于是本节提出了时空Hausdorff 距离来解决这一问题。

2.2 多维特征融合的相似性度量

Hausdorff 距离是描述2 个集合之间相似性的度量方法，可以用来度量轨迹间的相似性。优点在于不要求轨迹的采样率和点数量一致，可以度量整条轨迹的全局特征，当轨迹不完整时依然可以很好地进行相似性度量，具有较强的鲁棒性［25］。对于时空轨迹tri和trj，Hausdorff 距离的意义为tri中的位置信息点到trj中的位置信息点的最近距离中的最大值，定义tri到trj的定向Hausdorff 距离为

式中：{dist(li，k，lj，k)}表示2点之间的欧氏距离。以二维空间为例，

轨迹tri和trj的定向Hausdorff 距离h(tri，trj)如图2 所示。

图2 定向Hausdorff 距离示意图Fig.2 Schematic diagram of Hausdorff distance

通常情况下，采用双向Hausdorff 距离来定义相似度：

针对Hausdorff 距离只能度量轨迹位置特征而忽略速度、航向特征的问题，Pan 等提出了多因素Hausdorff距离［18］（Multi-Factor Hausdorff Distance，MFHD），同时考虑了轨迹的位置、速度以及航向信息，用多因素距离mfdist(pi，k，pj，k)代替欧氏距离来度量轨迹的多维特征相似度。

式中：wd、wv、wθ分别为位置、速度和航向信息的权重因子，权重因子的选择取决于应用场景，满足条件

然而，由于距离、速度和航向是量纲完全不同的信息，因此MFHD 很难确定一组行之有效的参数，存在一定的局限性。

事实上，目标在运动过程中的位置信息是目标的速度、航向等特征在时间上积累的结果，仅利用时间和位置信息即可很好地反映出目标在速度、航向等多维特征上的差异。因此考虑引入时空轨迹中的时间信息，通过构造时间滑窗，定义时空Hausdorff 距离（Spatio-Temporal Hausdorff Distance， STHD），由局部到整体逐步对时空轨迹进行相似性度量，利用时间积累后所有滑窗中最大的Hausdorff 距离来度量多维特征差异。时空Hausdorff 距离主要通过以下3 个步骤来构造：

步骤1确定时间滑窗长度τ。τ由参数ξ决定，ξ表示期望每个时间窗中信息数据点的个数。定义时间滑窗长度

步骤2根据滑动窗口长度τ对轨迹进行分段。若最后一段轨迹的时间不足一个滑窗长度τ，依然按照一个滑窗长度τ分段。分段后的时空轨迹tri和trj表示为

式中：tri，k和trj，k分别表示时空轨迹tri和trj中第k个时间窗中的轨迹段；M和N分别表示时空轨迹tri和trj中轨迹段的总数。

按时间顺序计算时空轨迹中2 个对应时间窗口中轨迹段的双向Hausdorff 距离H(tri，k，trj，k)，直到其中一个轨迹的时间滑窗终止。

步骤3计算时空轨迹tri和trj之间的时空Hausdorff 距离

这样就可以仅通过一个距离实现时空轨迹位置、航向以及速度等多维特征融合的相似性度量，同时也解决了多源时空轨迹数据采样率不稳定、数据多源异质等问题。

通过图3 中可以直观地看出，由于时间信息的引入，相比于Hausdorff 距离，对于位置相近而航向和速度不同的轨迹具有很好地区分度。

图3 时空Hausdorff 距离示意图Fig.3 Schematic diagram of STHD

3 多维特征融合的时空轨迹聚类

3.1 轨迹聚类

轨迹聚类是指通过最大化或最小化类内的相似性原则把轨迹数据集合分组为多个类或簇，同时用显式或隐式方法描述不同类或簇的过程，是轨迹数据挖掘的重要手段，对于揭示运动目标行为模式、侦查预警监视以及区域态势感知具有重要意义。目前轨迹聚类采用的主流算法包括K-means［9］、K-Medoids［10］、基于密度的带噪声的聚类方法［11］（Density-Based Spatial Clustering of Applications with Noise， DBSCAN）、密度峰值聚类［12］（Density Peaks Clustering， DPC）、谱聚类［13］（Spectral Clustering）等以及其改进算法。K-means 和K-medoids 算法简单易于实现，但是需要提前设置聚类参数，且迭代运算耗时过长。Wen 等［25］采用非参数统计方法核密度估计（Kernal Density Estimation， KDE），通过选取合适的核函数和滑动窗口自适应地解决轨迹聚类中存在的未知情况，但是选取核函数往往需要借鉴经验和专家系统。李旭东和成烽［26］基于密度峰值聚类，将密度峰值聚类中心作为经典轨迹，使挖掘出轨迹行为真实自然，但是该方法中截断距离的选择较为困难。Pan等［18］基于密度聚类的思想，利用MFHD 作为轨迹间的多维度相似性度量，提出了一种多维航迹聚类算法MTCA，但是该算法需要设置邻域半径ε和最小近邻航迹数量MinTRs，对于密度分布均匀的轨迹数据具有良好的效果，但是在目标行为模式分布不均匀的复杂数据集中存在虚警和漏警的问题。

3.2 时空轨迹k 最近邻决策聚类

密度峰值聚类［12］是近几年比较热门的聚类算法，该算法能够通过快速地搜索和发现密度峰值从而自动地发现簇中心，实现任意形状数据的高效聚类，在轨迹数据挖掘中得到了应用，但算法中截断距离的设定比较困难，并不能很好地适用于不同场景中。针对这一问题，我利用k最近邻思想，与密度峰值聚类中的决策图相结合，通过正整数参数k确定k最近邻距离，避免了截断距离参数的设定，同时采用时空Hausdorff 距离作为相似性度量，提出了多维特征融合的时空轨迹聚类算法（Spatio-Temporal Trajectory Clustering Algorithm， STTCA）。

首先介绍STTCA 的相关定义：

定义1在时空轨迹数据库TR 中，对于某一条时空轨迹tri，与之时空Hausdorff 距离第k近的时空轨迹定义为tri的k最近邻时空轨迹（k-Nearest Spatio-Temporal Trajectory），记为NSTT(tri，k)。

定义2时空轨迹tri与它的k最近邻时空轨迹的时空Hausdorff 距离定义为k最近时空Hausdorff 距离（k-Nearest Spatio-Temporal Hausdorff distance），记为

定义3定义最小k最近邻时空轨迹Hausdorff 距离

即在k最近时空Hausdorff 距离小于tri的所有时空轨迹中，与tri时空Hausdorff 距离最小的距离。特殊地，对于k最近时空Hausdorff 距离最小的那条时空轨迹，设置δ=max(STHD(tri，trj))。因此，只有那些k最近时空Hausdorff 距离小的轨迹，其最小k最近邻时空轨迹Hausdorff 距离才会远超正常值。

定义4以k最近时空Hausdorff 距离λi为横轴，最小k最近邻时空轨迹Hausdorff 距离δi为纵轴绘制出的散点图称为决策图。

定义5将最小k最近邻时空轨迹Hausdorff距离δi与k最近时空Hausdorff 距离λi的比值定义为决策值

式中：γi是进行聚类中心决策的重要参数。

定义6对于λi很大，即位于决策图靠右侧的轨迹，通过阈值THabnormal过滤，超过阈值的轨迹定义为异常轨迹。

定义7对于γi很大的轨迹，通过阈值THcenter决策，超过阈值的轨迹定义为中心轨迹。

定义8聚类完成后，对于某一时空轨迹tri，若它本身属于某一簇，但在其k最近时空Hausdorff 距离范围内存在属于其他簇的轨迹，那么定义其为可疑轨迹，属于异常轨迹的一种。

STTCA 算法流程如图4 所示，主要分为5 个步骤：

图4 STTCA 算法流程图Fig.4 Flow chart of STTCA algorithm

步骤1（轨迹相似性度量）计算时空轨迹数据库TR 中所有轨迹之间时空Hausdorff 距离，得到时空轨迹相似度矩阵。

步骤2（异常轨迹剔除）对时空轨迹数据库TR 中所有轨迹tri计算k最近时空Hausdorff 距离λi和与最小k最近邻时空轨迹Hausdorff 距离δi，以λi为横轴，δi为纵轴绘制出决策图，通过观察决策图中特征点的分布，在特征点密度很小的位置设定阈值THabnormal，λi大于该阈值的轨迹判定为异常轨迹，进行剔除。

步骤3（确定中心轨迹）计算决策值γi并按照降序排列，观察中心轨迹过渡到非中心轨迹的跳跃点，设定阈值THcenter，决策值大于阈值可判定为中心轨迹，从而确定轨迹簇个数。

步骤4（轨迹聚类）将剔除异常轨迹后剩余的轨迹进行分配。分配时，按照k最近邻时空轨迹Hausdorff 距离λi从小到大的顺序进行，将每条轨迹分配到k最近时空Hausdorff 距离比它小的轨迹中且与它时空Hausdorff 距离最近的轨迹所在的簇。

步骤5（可疑轨迹检测）对分配后的结果进行可疑轨迹检测，完成最终的目标行为规律挖掘。

4 仿真数据实验分析

本实验采用一个仿真的时空轨迹数据集模拟某预警监视区域中未知飞行器的飞行行为，利用本文所提算法和其他对比算法进行目标行为规律挖掘，并对实验结果进行对比分析。

4.1 仿真轨迹数据集

本实验使用Piciarelli 等［27］公开发表的目标航迹生成程序，该程序可以产生成簇的规律性轨迹以及随机的异常航迹，可设置轨迹簇密度分布参数，轨迹簇的个数，每一簇的轨迹条数，每条轨迹的信息数据点数，异常轨迹数等参数。用该程序产生二维空间中的目标轨迹数据，人为添加信息数据点中目标的时间、速度以及航向信息。本数据集共有1 500 条时空轨迹，每条时空轨迹包含20～60 个不等的信息数据点。在预警监视区域中总共设置了目标的10 种行为模式，即有10 个时空轨迹簇。

数据集的主要分布情况为主要为：在位置1附近生成4 簇轨迹，轨迹密度分布参数为0.5，每个轨迹簇有200 条轨迹，包括一簇标准行为、一簇2 倍速行为、一簇3 倍速行为和一簇反航向行为，用来验证算法对轨迹多维特征信息的挖掘能力；在位置2 附近生成2 簇轨迹，一簇为标准行为，密度分布参数为0.5，轨迹数量为200 条，另一簇为2 倍速行为，密度分布参数为0.6，轨迹数量为50 条，为了验证算法在高密度行为中挖掘低密度行为的能力；在位置3 附近生成2 簇轨迹，一簇为标准行为，密度分布参数为0.5，轨迹数量为100 条，另一簇为反航向行为，密度分布参数为0.6，轨迹数量为30 条，为了验证算法在高密度行为中挖掘更低密度行为的能力；在位置4 附近生成1 簇轨迹，密度分布参数为0.7，轨迹数量为100 条，为了检验算法对轨迹行为密度分布的鲁棒性；在位置5 附近生成1 簇轨迹，密度分布参数为0.7，轨迹数量为20 条，为了检验算法挖掘潜在行为的能力；同时还包括200 条随机生成的无规律异常轨迹，整个数据集如图5 所示。

图5 仿真轨迹数据集示意图Fig.5 Plot of all trajectories in dataset

4.2 算法复杂度

在STHD 算法中，对于度量含有n个信息数据点、N个轨迹段和含有m个信息数据点、M个轨迹段的2 条时空轨迹，本文算法复杂度为而MFHD 算法的复杂度为O(mn)，优于MFHD 算法。特殊地，当时间窗足够小时，STHD 近似于欧式距离度量，复杂度为O(n)，此时算法复杂度最小，对于速度和航向的度量效果最佳，但对轨迹整体的位置和形状度量效果最差；当时间窗足够大时，STHD 相当传统的Hausdorff 距离，复杂度为O(mn)，此时算法复杂度最大，对轨迹整体的位置和形状度量效果最佳，但无法度量轨迹的速度和航向信息。

在STTCA 算法中，对于含有n条时空轨迹的数据集，本文算法的复杂度为O(n2)，与MTCA 相当。

4.3 STTCA 算法实验结果

利用本文所提算法对仿真飞行器轨迹进行行为规律挖掘。采用时空Hausdorff 距离作为相似性度量，对于参数ξ的选择，考虑到STHD 的度量效果和算法复杂度，相似度计算消耗的时间随着ξ减小而减小，因此应当选择一个较小的值，但同时为了保证算法对轨迹整体的位置和形状度量效果，ξ又不宜过小，综合考量，设置ξ= 5 进行实验。采用STTCA 算法进行轨迹聚类，对于参数k的选择，考虑到算法复杂度，k越大，算法搜索近邻轨迹所消耗的时间越长，且参数k有一个重要意义：轨迹成簇和异常行为的分界点，即在某一位置聚集的轨迹数量超过k时就可以从海量轨迹数据中挖掘出该行为，聚集的轨迹数量未超过k时就不认为其为一种频繁出现行为，判定为无规律行为。因此为了挖掘出目标更多的行为规律，k的设定一般不宜过大。对于本节使用的仿真数据集，存在一簇包含20 条时空轨迹的行为规律和一簇包含30 条时空轨迹的行为规律，因此为了挖掘出这种低密度的行为，k的设置应当不超过20。分别设置k= 5， 15， 25， 35 进行实验，得到4 个λ-δ决策图如图6 所示。可以发现，当k= 5，15 时，都可以较好地挖掘出10 种行为规律的中心轨迹；而当k= 25 时仅可以发现9 种行为；当k= 35 时更是只能发现8 种行为。因此在后续的实验中，设置k最近邻参数k= 5。通过观察λ-δ决策图设置阈值THabnormal= 0.6，剔除异常轨迹。通过观察决策值降序图，如图7 所示，取决策阈值THcenter= 5，图中紫色虚线所示，提取出10 条中心轨迹。

图6 λ-δ 决策图Fig.6 λ-δ decision graph

图7 决策值示意图Fig.7 Plot of decision value

经过聚类和可疑轨迹判定，最终挖掘出飞行器运动的10 种行为规律，红色点代表目标轨迹起点，黑色点代表终点，相邻轨迹点间隔的长短表示目标的速度大小，如图8 所示。直接观察图5，很难发现目标轨迹的10 种行为规律，但是通过聚类结果图8，可以清晰直观地感知目标的10 种行为规律。

图8 挖掘出的10 种行为规律的中心轨迹Fig.8 Plot of central trajectories of excavated 10 behavior rules

为了证明本文算法的优越性，采用轮廓系数（Silhouette Coefficient）、精确度（Accuracy）和算法运行时间对算法进行评估，与轨迹相似性度量方法Hausdorff、MFHD 以及轨迹聚类算法MTCA、K-medoids、DPC 进行对比，采用不同参数进行多次重复实验，以多次试验结果的最优值作为最终结果，如表1 所示。

轮廓系数可以反映聚类出的簇的紧凑程度，轮廓系数越大表示同一簇中的目标行为越相近，不同簇中的目标行为越不相近，是一种内在的评估方法；精确率是指在同一个聚类结果簇中，与仿真数据集中基准确定的轨迹簇一致的轨迹占该聚类结果簇中所有轨迹的比例，是一种外在的评估方法。通过表1 可以看出本文算法与MTCA、K-medoids、DPC 挖掘出的行为模式在紧凑程度上相差不大。从算法实验结果精确度上来看，本文所提算法STHD 和STTCA 的组合性能最优，MTCA 算法由于对密度不均匀的行为挖掘能力较差，尤其是密度分布较低，轨迹数量较小的行为往往难以发现并被认定为异常轨迹；K-medoids 算法由于聚类中心随机且无法检测异常轨迹，会将不规律轨迹判定为某一种规律性行为，导致准确率不高；Hausdorff 只能度量轨迹之间的位置信息而无法度量速度和航向信息，因此造成整体的挖掘算法准确率不高。从算法运行时间上看，除了K-medoids 聚类迭代消耗大量时间外，其他的算法时间消耗主要集中在轨迹相似性度量上，STHD 在时间复杂度上低于MFHD和Hausdorff，避免了不必要的信息点距离运算，整个算法的运行时间大大降低。

表1 算法性能比较Table 1 Comparison of algorithm performance

5 实测数据实验验证

本实验选取了2006年1月某天旧金山机场雷达检测出的所有航班飞行轨迹，经过初步筛选，剔除部分无效轨迹，总共有1 017 次飞行时空轨迹，每条时空轨迹包含x，y和z坐标，其中z为飞行高度，以及雷达信息数据点的时间信息等，所有轨迹如图9 所示。

利用本文所提算法对航班飞行轨迹进行行为规律挖掘。采用时空Hausdorff 距离作为相似性度量，设置参数ξ= 5；采用STTCA 算法进行轨迹聚类，设置k最近邻参数k= 5 进行目标行为规律挖掘实验。通过观察λ-δ决策图设置阈值THabnormal=1×104m，剔除异常轨迹，如图10 所示，红色圈内的轨迹被认定为无规律轨迹。

图10 实测数据λ-δ 决策图Fig.10 λ-δ decision graph of measured data

通过观察决策值降序图，如图11 所示，取决策阈值THcenter= 10，图中紫色虚线所示，提取出12 条中心轨迹。经过聚类和可疑轨迹判定，最终挖掘出航班飞行的12 种行为规律，如图12 所示，其中红色点表示目标行为的起点，黑色点表示目标行为的终点。实测数据实验表明，本文算法在真实应用场景中具有良好的规律挖掘性能，在智能监视任务中具有较好的应用前景。