随着国内外经济的不断发展，船舶数量不断增加，并呈现大型化、高速化和智能化的发展趋势^[1]。沿海港口、河口段等水域船舶通航密度不断提高，水域通航环境更加复杂，对水域通航管理能力提出了更高的需求。船舶自动识别系统(automatic identification system, AIS)作为船舶运动信息的载体为相关研究提供了海量的数据。对AIS数据进行挖掘，聚类分析船舶轨迹信息，提取水域船舶主要航路，可为海事监管及船舶交通管理系统(vessel traffic service，VTS)智能化提供支持^[2]。

在水上交通方面，国内外学者在轨迹聚类开展了研究。肖潇^[3]在AIS数据的基础上，通过选取航迹特征点，对船舶轨迹进行分段，利用基于密度的聚类算法(density-based spatial clustering of applications with noise，DBSCAN)对船舶轨迹进行聚类，得到了水域船舶主要航路。马文耀等^[4]提出基于单向距离的轨迹相似度度量方法，采用谱聚类对水域船舶轨迹进行聚类分析，得到水域船舶行为模式。这些算法提高了对数据中噪声和孤立点的容忍度，但由于DBSCAN等聚类算法没有考虑轨迹整体形状特征，在聚类过程中对轨迹进行了划分，难以有效地提取整段轨迹。当数据分布密度变化较大时，该类聚类方法的结果往往较差，也有学者针对整段轨迹展开聚类研究^[5-6]，但算法对计算机内存等资源消耗较大。为了克服以上算法中的问题，本文提出了一种结合Douglas-Peucker算法与谱聚类的船舶轨迹快速自适应谱聚类算法。

1 船舶轨迹数据预处理 1.1 AIS数据解码与清洗

AIS数据不仅记录船型、船舶尺度等静态信息，而且以每隔2 s~6 min的时间间隔发射船位、航速和航向等动态信息，相关数据的编码方式遵循ITU-RM.1371-4和IEC61162协议分为明码和暗码报文。本研究编写AIS解析程序，从原始数据中提取研究所需的船舶静态数据(船舶水上移动业务标识、船舶尺度等)、动态数据(位置、速度、航向等)。此基础上，过滤异常船速、航向点，去除错误、重复数据。实际AIS数据中，存在一艘船舶轨迹数据包含两艘并不连续的轨迹的情况，为了避免对聚类结果造成干扰，需要对这些轨迹进行分割处理。

1.2 基于Douglas-Peucker(DP)算法的船舶轨迹数据压缩

Douglas-Peucker算法最早由David Douglas和Thomas Peucker^[7]于1973年提出，是经典的线状要素简化算法。算法思路为：对一条由点构成的曲线，连接曲线首尾点得到一条新的直线；计算曲线上所有点到该直线的距离并求取最大值；比较距离最大值d_max与算法阈值d_th:若d_max＜d_th，则删去该曲线除首尾点以外所有的点; 若d_max≥d_th, 则保留d_max对应的点，并根据该点将曲线分成两个部分；对两部分曲线分别重复上述过程，直至曲线无法进行分割，所得即为简化后的曲线。

DP算法的阈值设置方面，以某船舶水上移动业务标识为"209XXXXXX"的轨迹为例，如图 1所示，较小的阈值即可取得较高的压缩比。阈值为20 m时，轨迹压缩比为98.13%。随着阈值增大，压缩比无显著增加。图 2给出某船舶轨迹经阈值在100 m的DP压缩处理压缩前、后的对比。

船舶原始轨迹数据中包含1 545个数据点，压缩后保留18个数据点，压缩比为98.83%。从图 2可以看出，压缩后的数据仍反映了该轨迹的整体形状特征。为了让压缩后的轨迹尽量保留原轨迹的形状特征，阈值的设置以不超过区域最窄航道宽度的四分之一为宜。

张树凯等^[8]应用DP算法，提高了电子海图显示与信息系统(electronical chart display and information system, ECDIS)船舶轨迹显示效率，并验证了该算法能够在保留轨迹特征点的情况下以较低的失真度对轨迹进行压缩。因此本文采用该算法对船舶轨迹数据进行压缩。

2 相似度度量函数 2.1 Hausdorff距离

当前基于距离的轨迹相似度度量方法多采用Euclidean距离、Hausdorff距离、结构相似度距离等^[3]。考虑空间目标的形状差异与相对位置差异，德国数学家Felix Hausdorff提出一种新的距离度量方法，即Hausdorff距离^[9]。假设两个点集A={a₁, a₂, a₃, …, a_i}与B={b₁, b₂, b₃, …, b_j}，则A与B的Hausdorff距离可由下式得到

$ {D_{{\rm{hausdorff}}}} = \max \left( {d\left( {A, B} \right), d\left( {B, A} \right)} \right) $

(1)

$ d\left( {A, B} \right) = \mathop {\max }\limits_{{a_i} \in A} \left( {\mathop {\min }\limits_{{b_j} \in B} \left\| {{a_i} - {b_j}} \right\|} \right) $

(2)

$ d\left( {B, A} \right) = \mathop {\max }\limits_{{b_j} \in B} \left( {\mathop {\min }\limits_{{a_i} \in A} \left\| {{b_j} - {a_i}} \right\|} \right) $

(3)

式中：d(A, B)、d(B, A)分别为A到B、B到A的单向Hausdorff距离，该值可通过穷举计算点集间最小距离中的最大值得到; ‖‖表示点之间的距离。船舶位置数据是以经纬度的形式通过AIS进行播发，本研究中采用球面距离以提高计算精度。

相对于传统的距离度量，Hausdorff距离未对数据集进行插值，进而避免了对轨迹数据添加噪声。另一方面，该方法可以度量不同长度船舶轨迹之间的hausdorff距离，对由于播发时间间隔不一致以及长度不同的船舶AIS轨迹数据具有良好的适应性。因此，本研究采用Hausdorff距离进行相似度函数的构造。

2.2 自适应尺度参数的相似度度量函数

一般情况下，轨迹之间的相似度函数如下

$ s\left( {A, B} \right) = {{\rm{e}}^{ - \frac{{d{{\left( {A, B} \right)}^2}}}{{2{\sigma ^2}}}}} $

(4)

式中：s表示轨迹间的相似度；d(A, B)表示轨迹间的距离，这里为Hausdorff距离；σ为尺度参数，表示相似度随着距离的增大而衰减的程度^[4]。

通常尺度参数的选取非常繁杂，需要对参数进行不断尝试才能取得满意的结果。为了改进这一缺陷，一些学者针对尺度参数的选取展开了研究：卜德云^[10]基于类似核选取的方法提出了自动选取尺度参数σ的方法，但该算法需要进行大量重复计算，计算开销较大；基于"local scale"^[11]的思想，Chen等^[12]在设计分布式谱聚类算法的过程中，提出可以用轨迹相对于其他所有轨迹距离的平均值或中值代替尺度参数σ。本文根据文献[12]，采用下式所示的相似度度量函数：

$ s\left( {A, B} \right) = {{\rm{e}}^{ - \left[{\frac{{d{{\left( {A, B} \right)}^2}}}{{2{\sigma _{{A^\sigma }B}}}}} \right]}} $

式中：σ_A、σ_B分别为轨迹A、B与其他轨迹间Hausdorff距离的平均值，用以代替传统需要人工确定的尺度参数。

3 船舶轨迹谱聚类

谱聚类(spectral clustering, SC)源于谱图划分理论^[13]，该算法通过对聚类对象的相似度矩阵进行特征分解，求解矩阵的特征值或特征向量，选择合适的特征向量对数据进行聚类。采用谱聚类的优点在于可在任意形状的样本空间上聚类, 具有全局最优性^[10]，并且便于计算机实现。算法主要步骤如下：

1) 对n个聚类对象，构造n阶相似度矩阵A^n×n：

$ {\mathit{\boldsymbol{A}}_{i, j}} = \left\{ \begin{array}{l} 0, \;\;\;\;\;\;0 < i = j \le n\\ {S_{i, j}}, \;\;\;\;0 < i \ne j \le n \end{array} \right. $

式中S_{i, j}为轨迹i与j的相似度。

2) 构造对角矩阵D^n×n，对角线上元素为相似度矩阵某行或列元素的和：

$ {\mathit{\boldsymbol{D}}_{i, i}} = \sum\limits_j^n {{\mathit{\boldsymbol{A}}_{i, j}}} $

3) 计算并归一化拉普拉斯矩阵L^n×n：

$ \begin{array}{l} \;\;\;\;\;{\mathit{\boldsymbol{L}}^{n \times n}} = {\mathit{\boldsymbol{D}}^{n \times n}} - {\mathit{\boldsymbol{A}}^{n \times n}}\\ {\mathit{\boldsymbol{L}}^{n \times n}} = {\mathit{\boldsymbol{D}}^{n \times n - \frac{1}{2}}} \cdot {\mathit{\boldsymbol{L}}^{n \times n}} \cdot {\mathit{\boldsymbol{D}}^{n \times n - \frac{1}{2}}} \end{array} $

4) 对得到的拉普拉斯矩阵L^n×n进行特征分解，计算特征值与特征向量，将所得矩阵K个特征值从小到大排列，并将对应特征向量构成一个新的N×K阶矩阵，采用K-Means算法对特征向量进行聚类，得到轨迹聚类结果。

本研究目的在于提出一种能够对船舶轨迹数据进行快速聚类分析，提取水域船舶主要航路的聚类方法，算法逻辑流程图如图 3所示。

4 实验与验证

长江口作为整个长江水域的"咽喉"部位，主要由北槽(深水)航道、南槽航道、北港航道以及北支航道组成，连接我国长江内河航运与远洋运输，船舶交通非常繁忙，著有"三级分汊、四口入海"之称。本文基于C#语言与visual studio平台，编写船舶轨迹自适应谱聚类算法的程序并进行了实验。实验采用2010年3月16日长江口水域(30°58'11″N~31°18'56″N，121°43'19"E~122°32'24"E)实际AIS数据为测试样本。长江口船舶轨迹如图 4所示。

AIS数据中包含表示船舶类型的字段，通过数据分析发现，水域船舶中包含一定数量从事拖带、引航、水下作业等活动的特殊船舶，这些船舶与从事客货运输船舶的活动特征不同，例如在小范围水域内往返多次等。为了避免对聚类造成影响，这些船舶的轨迹需要从样本中剔除。

经过对1 353 789个原始AIS数据进行解码、清洗、剔除、分割等操作，得到了359艘船舶轨迹。在DP压缩的阈值选取方面，根据上海海事局公布的上海港通航环境，该水域主要航路的宽度在500~2 315 m且双向通航，本研究选取100 m作为阈值，对轨迹应用FSSC算法进行聚类。经过两轮聚类，得到最终准确的轨迹分类，根据结果将AIS数据分组在谷歌地图上进行展示。

如图 5所示，经过对聚类数据的分析与整合，算法得到了长江口水域船舶主要的7条航路，分别为：航路1(图 5(a))经由长江口主航道与长江口1#灯船连接长江与北部海域；航路2(图 5(b))经由长江口主航道与长江口1#灯船连接长江与东部海域；航路3(图 5(c))经由长江口主航道与长江口1#灯船连接长江与南部海域；航路4(图 5(d))经由南槽航道与长江口2#灯船连接长江与北部海域；航路5(图 5(e))经由南槽航道与长江口2#灯船连接长江与东部海域；航路6(图 5(f))经由南槽航道上段与南支航道连接长江与南部海域；航路7(图 5(g))经南槽航道与长江口1#灯船连接东部与南部海域。

在首轮聚类中，得到了航路1、4、6，而航路2、3、5、7需要通过第2轮FSSC聚类得到。这是由于航路2与3、5与7包含的轨迹间存在较长的相似航迹段，差异较小，在总体样本中不足以被准确划分。在第一轮聚类的基础上，对没有准确聚类的数据应用二次FSSC后得到了准确结果。实验对比了同样采用距离度量轨迹相似性的层次聚类法与谱聚类在相同DP压缩阈值和初始聚类数下第一轮聚类的结果。谱聚类算法成功聚类得到3条航路，层次聚类法仅得到1条。从聚类效率上看谱聚类算法优于层次聚类。

根据2008年颁布的《长江口船舶定线制》，长江口船舶定线制的主要组成部分为2个警戒区、5条分隔带及5条通航分道，并与其相关的锚地和引航作业点共同构成完整的船舶交通体系，如图 6所示。由该定线制规定与南支航道一同将长江口水域划分为6条主要航路。由此可知本文提出的船舶轨迹聚类算法能够较为准确地对船舶轨迹进行聚类，并且结果符合水域实际交通状况。算法的计算速度方面，采用DP算法后，计算时间从原本的数小时下降到不足20 min，算法的计算速度有明显提高。

5 结论

1) 该算法能够对不同航路进行准确区分并聚类；

2) 算法计算速度明显快于普通聚类方法。

在具体细节上，该算法难以一次得到准确的聚类结果，下一步研究将针对这一问题将开展自适应多轮FSSC聚类的研究。DP算法应用中的阈值设定仍需要进行进一步完善，算法的计算速度可通过采用并行计算等技术进一步提高。