0 引　言

随着航海数字化的发展，船舶自动识别系统（Automatic Identification System，AIS）得到了广泛应用。AIS数据中蕴涵了大量船舶运动特性^[1]。因此，AIS数据常被用于交通流特征提取^[2]、航路规划^[3]、轨迹预测^[4]、异常行为检测^[5]等研究，这些研究离不开对船舶行为模式分析和交通流分布特征的了解，轨迹聚类是实现这类研究的有效手段。

目前常用的船舶轨迹聚类方法有K-Medoids聚类^{[6 – 7]}、基于密度的空间聚类（Density-Based Spatial Clustering of Applications with Noise，DBSCAN）聚类^{[8 – 9]}、谱聚类^[10]和基于密度的轨迹聚类^[11]（Density-Based Trajectory Clustering of Applications with Noise，DBTCAN）等算法。K-Medoids聚类的关键是设置合适K值，但K的选择易受噪声的影响^[12]。DBSCAN对噪声具有较好的鲁棒性，但2个关键参数需要根据经验选择，吸引了较多的学者深入研究^[13]。谱聚类可以处理稀疏数据，但参数较多，性能也依赖于相似矩阵^[10]。上述3种聚类方法忽略了位置点间的时空关联性，DBTCAN以整条轨迹作为研究对象可以更完整的提取和分析船舶运动模式。

度量点和轨迹之间是否相似是轨迹聚类的关键。常用的度量方法是使用距离进行表征，如欧式距离^[14]（Euclidean Distance，ED）、豪斯多夫距离^{[11, 15]}（Hausdorff Distance，HD）、弗雷歇距离^[16]（Fréchet Distance，FD）和动态时间规整距离^[17]（Dynamic Time Wraping，DTW）。ED多适用于轨迹点。HD可以较好的适用于长度和轨迹点数量不同的轨迹之间，但易受轨迹中存在异常离群轨迹点的影响。如果轨迹中存在异常离群轨迹点，就会引起非常大的距离变化，对噪声非常敏感且无法区分轨迹方向。FD和HD都是基于形状的距离度量方式，因此同样存在对轨迹中噪声敏感的缺点。DTW在噪声处理方面优于HD，但DTW计算复杂度高、效率低。BAI等^[17]提出Fast-DTW以提高效率，但仅利用了位置信息，无法区分轨迹方向。多维DTW考虑了空间和运动特征差异^[18]，但在计算DTW过程中每个轨迹点都会多次被重复使用，所以在处理极端匹配的情况（如多对一和一对多）时会导致距离过度失真，使得轨迹间的距离发生实质性变化。

通过上述分析发现，现有的轨迹聚类方法普遍存在确定最佳参数困难的问题；常用的通过距离表征轨迹相似的度量方法对轨迹中的噪声敏感，无法给出相似轨迹之间的相似程度和区分轨迹的方向。针对以上问题，首先，融合空间位置和航向特征构建航向约束的分段路径距离（Course over Ground Segment-Path Distance，CSPD），并定义轨迹相似度，能连续地描述和表征轨迹的相似程度，受轨迹中异常噪声点的影响较小，且可以区分不同方向的轨迹；其次，根据轨迹相似度分布特征和聚类评价指标构建自适应确定参数流程，参考文献[11]的DBTCAN方法，本文提出一种基于相似度和密度的抗噪声轨迹聚类（Degree of Similarity and Density Based Trajectory Clustering of Applications with Noise，DS-DBTCAN）方法。

1 轨迹相似度度量方法 1.1 轨迹空间距离

空间距离是度量轨迹相似度的主要依据。常用的距离存在轨迹空间距离的度量结果易受轨迹长度差异和异常噪声数据等问题的影响。因此，提出一种改进的轨迹空间距离，适用不同长度的轨迹且受轨迹中偶然变化的影响较小。给定轨迹T₁={p₁,…,p_i,…,p_m}和T₂={p₁,…,p_j,…,p_n}。

定义1　点到子轨迹段投影距离（Projected Distance，PD）是点到另一条轨迹的子轨迹段的最短欧式距离，即：

${D_{pd}}\left( {p,T} \right) = \mathop {\min}\limits_{{s_j} \in T} {d_{\min }}\left( {{p_i},{s_j}} \right)。$

(1)

式中：s_j为轨迹T的子轨迹段；d_min(p_i,s_j)为点到轨迹段的最短欧式距离，如图1(a)所示。

定义2　分段路径距离（Segment-Path Distance，SPD）是指一条轨迹上各轨迹点到另一条轨迹投影距离的平均值，即：

${D_{spd}}\left( {{T_1},{T_2}} \right) = \dfrac{1}{m}\sum\limits_{{p_i} \in {T_1}} {{D_{pd}}\left( {{p_i},{T_2}} \right)}，$

(2)

${D_{spd}}\left( {{T_2},{T_1}} \right) = \dfrac{1}{n}\sum\limits_{{p_j} \in {T_2}} {{D_{pd}}\left( {{p_j},{T_1}} \right)}。$

(3)

式中：m、n分别为轨迹T₁和T₂中轨迹点数量；D_spd(T₁,T₂)为T₁中所有点到T₂的投影距离的平均值；图1(b)为SPD计算示例。

SPD是各轨迹点到另一轨迹的投影距离的平均值来度量轨迹间距离，削弱了噪声点的影响。但SPD是非对称距离，故本文将两轨迹之间SPD距离的最小值作为轨迹空间距离，即：

${D_d}\left( {{T_1},{T_2}} \right) = \min \{ {D_{spd}}({T_1},{T_2}),{D_{spd}}({T_2},{T_1})\}。$

(4)

1.2 方向差异性

在轨迹空间距离的基础上，T₁中每个轨迹点p_i都可以在T₂上找到与之相对应的投影点q_i，计算轨迹点与其对应投影点之间运动方向的差异性。对于轨迹T={(m₁,t₁,x₁,y₁,c₁),…,(m_i,t_i,x_i,y_i,c_i),…,(m_n,t_n,x_n,y_n,c_n)}，相邻两轨迹点间的轨迹段矢量方向计算如下式：

${ {C_{{P_i}{P_{i + 1}}}} = \left\{ \begin{gathered} {\sin ^{ - 1}}\Biggr( \left( {{x_{i + 1}} - {x_i}} \right)/ \\ {\sqrt {{{\left( {{x_{i + 1}} - {x_i}} \right)}^2} + {{\left( {{y_{i + 1}} - {y_i}} \right)}^2}} } \Biggr),{x_1} \leqslant {x_2}，\\ 360^\circ - {\sin ^{ - 1}}\Biggr( \left( {{x_{i + 1}} - {x_i}} \right)/ \\ {\sqrt {{{\left( {{x_{i + 1}} - {x_i}} \right)}^2} + {{\left( {{y_{i + 1}} - {y_i}} \right)}^2}} } \Biggr),{x_1} > {x_2}。\\ \end{gathered} \right. ,i = 1,2,...,n。}$

(5)

轨迹之间方向差异性的计算具体描述如下：

1）给定2条轨迹T₁和T₂，p_i为T₁上的轨迹点，q_i为轨迹点p_i在T₂上的投影点，如果q_i是T₂上的轨迹点，则运动方向为该点航向Cq_i，如果q_i不是T₂上的轨迹点，则以q_i所在轨迹段依据式(5)计算的矢量方向Cp_i p_i+1表示q_i的方向，方向差异性计算如图2所示。

2）计算T₁上的各轨迹点p_i与其在T₂上的投影点q_i运动方向的差异性d_PC(p_i,T₂)，即：

${ {d}_{PC}\left({p}_{i},{T}_{2}\right) = \left\{ \begin{array}{l}\dfrac{1-\mathrm{cos}\text{(}C{p}_{i}-C{q}_{i}\text{)}}{2}，{q}_{i}{\text{是}}T2{\text{上轨迹点}}，\\ \dfrac{1 - \mathrm{cos}\text{(}C{p}_{i} - C{p}_{i+\text{1}}{p}_{i}\text{)}}{2}，{q}_{i}{\text{不是}}T2{\text{上轨迹点}}。\end{array} \right.}$

(6)

3）轨迹T₁相对于T₂的方向差异性d_TC(T₁,T₂)，即：

${d_{TC}}\left( {{T_1},{T_2}} \right) = \frac{1}{n}\mathop \sum \limits_{{p_i} \in {T_1}} {d_{PC}}\left( {{p_i},{T_2}} \right)。$

(7)

4）轨迹T₁和T₂间运动方向的整体差异性D_C，D_C∈[0,1]，值越小则轨迹T₁和T₂方向越接近一致。

${D_C}\left( {{T_1},{T_2}} \right) = \max\left\{ {{d_{TC}}\left( {{T_1},{T_2}} \right),{d_{TC}}\left( {{T_2},{T_1}} \right)} \right\}。$

(8)

1.3 轨迹相似度

空间位置和运动方向是度量轨迹相似度的2个重要特征。考虑到空间距离是度量轨迹相似度的主要因素，因此，将方向差异性作为轨迹空间距离的权重来构建CSPD距离，即：

$\text{CSPD}({T}_{1},{T}_{2})=\frac{1}{1-{D}_{C}({T}_{1},{T}_{2})} \cdot {D}_{d}({T}_{1},{T}_{2})。$

(9)

式中：1/(1-D_C(T₁,T₂))表示轨迹T₁和T₂之间的方向差异性对空间距离的修正权重，该权重随着D_C的增加（方向的差异性增大）而单调递减。若2条轨迹方向完全相同，即D_C =0，该权重为1；若2条轨迹运动方向完全相反，即D_C =1，该权重趋向于无穷大，能够最大程度上区分空间相近但方向不同的轨迹。

目前度量轨迹相似的方法是基于轨迹距离来描述的，并且需要有互相参照的对象进行对比，其结果不够简单连续直观和理解。例如，只能衡量2个轨迹和同一个轨迹之间，从而判断哪2个轨迹更相似，但轨迹距离单调递增，无法单纯的判断和表征2个轨迹之间的相似程度。因此，本文基于轨迹CSPD距离定义了轨迹相似度，首次连续直观地表征和描述两轨迹之间的相似程度，易于理解，且结果更方便进行归一化处理。船舶轨迹相似度函数公式如下：

$S({T}_{1},{T}_{2})=\frac{{k}^{2}}{{k}^{2}+ \text{CSPD}({T}_{1},{T}_{2})}。$

(10)

式中：k为常量参数，轨迹之间的相似度函数随着距离的变化曲线（k=5时）如图3所示。

本文分别采用HD、DTW、FD、改进的HD^[15]以及本文提出的基于CSPD距离的轨迹相似度方法，模拟设计了4种类型的轨迹进行相似度的度量实验。如图4所示，结果表明基于CSPD距离的轨迹相似度方法能更连续直观的描述和表征轨迹的相似程度，不仅受轨迹中噪声轨迹点的影响较小，具有较好的鲁棒性，还可以有效区分不同运动方向的轨迹。

2 基于相似度和密度的抗噪声轨迹聚类

首先，对轨迹数据进行清洗和预处理；然后，基于融合位置和航向特征的CSPD距离定义轨迹相似度；最后，根据聚类评价指标构建自适应参数确定流程，利用DS-DBTCAN算法对轨迹聚类分析并提取特征轨迹，轨迹聚类流程如图5所示。

2.1 DS-DBTCAN算法

DBTCAN算法可以直接对轨迹进行聚类。基于CSPD的船舶轨迹相似度方法是以相似度进行衡量描述和表征轨迹之间的相似程度，所以DS-DBTCAN的输入参数分别为相似度阈值S_min和相似轨迹最低数量L_min。DS-DBTCAN的基本原理是选取轨迹集中的任意一条轨迹，找到与该轨迹相似度大于等于S_min的所有轨迹：若轨迹数量小于L_min，则该轨迹被标记为噪声轨迹；若轨迹数量大于等于L_min，则该轨迹为核心轨迹。接着判断与这条轨迹相似度大于等于S_min的其他轨迹是否为核心轨迹，如果是，则这些轨迹属于一个簇C。再对其他未被选中的轨迹进行判断，直到所有轨迹判断完毕完成聚类为止。

2.2 自适应确定最佳聚类参数

实现参数的自适应确定，关键在于确定合适的候选阈值参数列表集合。首先基于轨迹集自身的相似度分布特性，利用K平均最近邻（K-average nearest neighbor，KANN）算法^[19]和数学期望法分别生成DS-DBTCAN算法的输入候选参数S_min和L_min列表集合。由聚类评价指标得分最高的参数S_min和L_min则为最佳的聚类参数，即该参数所得出的聚类簇为最优聚类簇。步骤如下：

步骤1　采用KANN算法生成候选相似度阈值S_min列表。KANN算法的基本思想是通过计算轨迹集中每条轨迹与其第K条最近邻轨迹之间的K最近相似度，并对所有轨迹的K最近相似度求均值得到轨迹集的K平均最近相似度。

1）计算轨迹集D的相似度矩阵，即：

${\boldsymbol{D}_{l \times l}} = \left\{ {H\left( {i,j} \right)|1 \leqslant i \leqslant l,1 \leqslant j \leqslant l} \right\}。$

(11)

式中：l为轨迹集D所包含的轨迹数量。

2）对相似度矩阵D_l×l中的每行元素降序排列，第K列元素构成所有轨迹的K最近相似度向量D_K。

3）对D_K中的所有元素求平均值，得到K平均最近邻相似度$\overline {{D_K}}$，将其作为候选的S_min。对所有的K值进行计算，则得到长度为n的候选S_min列表，表示为D_s：

${D_s} = \left\{ {\overline {{D_K}} \mid 1 \leqslant K \leqslant n} \right\}。$

(12)

步骤2　采用数学期望法生成候选L_min列表。计算所有轨迹的相似度大于等于S_min的轨迹数量的数学期望值，作为轨迹集D的长度为m的候选相似轨迹最低数量阈值L_min列表。

$L_{\min}=\dfrac{1}{n} \sum_{i=1}^{n} P_{i}。$

(13)

式中：P_i为与第i条轨迹相似度大于等于S_min的轨迹数量。

步骤3　确定最佳相似度阈值S_min和相似轨迹最低数量阈值L_min。依次将S_mini和L_minj输入DS-DBTCAN算法对轨迹集进行聚类分析，可以得到不同K值下的聚类簇数量，并计算每个聚类簇的综合聚类评价指标，选择指标得分最高的类簇作为最优聚类簇，自适应确定最佳参数流程如图6所示。

2.3 聚类评价指标

聚类评价指标用于评估聚类性能和结果，因此将轮廓系数（Silhouette Coefficient，SC）和邓恩指数（Dunn Validity Index，DVI）综合起来评价聚类的性能和效果。综合聚类评价指标（Comprehensive Clustering Metrics，CM）的公式如下：

${CM}=\frac{{SC}'+{DVI}'}{2} 。$

(14)

式中：${SC}'$、${DVI}'$分别为标准化后的SC和DVI。CM值越大，表明聚类效果越好。

2.4 特征轨迹提取

特征轨迹是同类簇内船舶运动行为的最佳表达形式，是每个轨迹聚类簇中具有代表簇内整体特征信息的轨迹，因此特征轨迹可以表示习惯航线和交通流信息。通过扫描线法^[20]对聚类后的船舶轨迹簇，提取具有代表性的轨迹来描述船舶的整体运动特征并可视化。首先建立一条垂直于轨迹簇的整体运动方向的扫描线，记录扫描线与轨迹的交点，如果交点数量超过阈值，计算平均坐标等信息形成特征轨迹中心线，如图7所示。

3 方法验证与实例分析 3.1 仿真实验结果与分析

为了验证DS-DBTCAN算法的有效性，实验模拟了不同运动方向的轨迹进行聚类，分别以DTW、改进的HD^[15]和基于CSPD距离的轨迹相似度方法，计算基于不同轨迹相似性度量的方法得到的聚类结果的F值^[21]。结果如表1所示，采用基于CSPD距离的轨迹相似度的方法时，聚类结果的总体精度优于其他的度量方法。K=32时综合聚类指标得分最高，此时所对应的相似度阈值S_min=0.82、相似轨迹最低数量阈值L_min=10，能获得9个类簇，取得较为理想的聚类结果。如图8所示，可以发现使用DTW和改进的HD方法时将几何形态比较相似和空间位置距离非常相近的轨迹错误地聚成一类，尤其是部分轨迹存在相互重叠，错误更明显。

3.2 真实数据实验结果与分析

长江口是我国航运密度最大的区域之一，选取长江口水域（30.85°～31.30°N，121.7°～122.5°E）2020年3月376艘船舶308838个轨迹点，轨迹分布如图9所示，进行轨迹聚类和特征轨迹提取，验证算法对船舶轨迹聚类的适应性。

利用DS-DBTCAN方法对长江口AIS数据进行分析，得到聚类簇数量与K值的关系，如图10(a)所示。从K=73开始，类簇数量保持不变且趋于稳定状态，直到K>86后，类簇数量逐渐变小直至趋于稳定不再变化，而此时对应的综合聚类评价指标CM先增大并且不断波动，最后趋于稳定。在K=83时，CM值最高，说明此时的聚类结果较优，对应的K平均最近邻相似度即为最佳的相似度阈值S_min=0.76。利用式(13)计算大于等于该相似度阈值S_min所对应的相似轨迹数量，即可得到最佳的相似轨迹最低数量阈值L_min=6，能获得较好的聚类效果，如图10(b)所示，聚类出8个轨迹类簇，同一航路的颜色不同的2个类簇表示轨迹的方向不同，代表了航路的上下行交通流。

根据图10(b)可以发现聚类结果成功识别出海图中所示航道的上行与下行航线，与水域实际航路分布相吻合，类簇1和类簇2是经由南槽航道上段与南槽支航道连接长江与南部海域的航路，类簇3和类簇5是洋山港方向的船舶穿越航道进出南槽航道的习惯航路，类簇4和类簇6为南槽航道连接长江与东部海域的航路，类簇7和类簇8为长江口北槽主航道的航路。

根据特征轨迹的定义和提取方法，提取出具有方向的8条特征轨迹，如图11所示，航道的上下行特征轨迹均能清晰提取出来。将特征轨迹与标准航线的相似度作为误差判断，提取电子海图标准航道的分隔线与其航道2条边界线的平均值，即可得出航道中心线作为标准航线。特征轨迹与标准航线的相似度为0.78、0.83、0.92、0.86，结合实际水域交通状况和相关航路信息，提取出的特征轨迹误差较小，效果比较理想和真实航道情况相吻合。

4 结　语

本文提出基于CSPD距离的轨迹相似度相比于传统距离判别轨迹相似的方法，不仅可以连续地表征轨迹之间的相似程度，而且能区分不同方向的轨迹，对轨迹中噪声点具有鲁棒性。模拟轨迹数据聚类实验表明，基于CSPD距离的轨迹相似度方法得到的F值精度优于使用其他度量方法的聚类结果。以长江口水域AIS数据为例，DS-DBTCAN算法能够在大量复杂船舶轨迹中对特征相似的轨迹聚类出多个不同向的轨迹簇，区分出的上下行航路与实际交通流情况基本吻合，可为相关部门在航路规划和开展定线制工作提供一定的参考。