基于大规模船舶轨迹数据的航道边界提取方法

徐垚 李卓然 孟金龙 赵利坡 温建新 王桂玲

摘 要：传统的道路数据获取方法成本高、更新慢等无法适用于海洋航道的获取，从众源轨迹数据中提取道路或航道信息具有成本低、更新快等特性，然而，由于船舶轨迹数据噪声多、数据量大、不同区域分布不均使得航道边界提取面临较大挑战。针对该问题，提出一种基于大规模船舶轨迹数据进行航道边界提取的方法。首先对大规模的船舶轨迹数据进行并行化去噪、插值、轨迹分段;然后，基于并行化及基于Geohash编码的空间聚类，将轨迹数据化简为多个方形区域的点集数据;其次，对其进行窗口划分，对传统的NiBlack方法进行扩展，提出SpatialNiBlack算法，对方形区域进行航道识别;最后，提出一种新的提取算法del-alpha-shape，基于航道识别结果获得航道边界。理论分析与实验结果表明，所提方法在最大密度值是200，最小密度值是10，窗口长和宽分别为5和5时，可同时达到86.7%的准确率和79.4%的召回率。实验结果表明，该方法可以从大规模的轨迹数据中提取有价值的航道边界，是一种有效的航道提取方法。

关键词：轨迹数据;自动识别系统;时空大数据;Delaunay三角网;航道提取

中图分类号： TP319

文献标志码：A

Abstract： The traditional road information extraction method is high-cost and slow-update. Compared with it， road or marine lane information extraction from crowdsourcing trajectory data is low-cost and easier to update. However， it is difficult to extract lane boundary due to vessel trajectory data with high noise， large data volume and uneven distribution across different regions. To solve this problem， an extraction method of marine lane boundary from exploiting trajectory big data was proposed. Firstly， the parallelized denoising， interpolation and trajectory segmentation for trajectory big data was conducted. Then， based on parallelization and Geohash-encoded spatial clustering， trajectory data was simplified into multiple square regions. The regions were divided and the NiBlack method was extended as SpatialNiBlack algorithm to recognize regions on lane. Finally， based on the filtering results， del-alpha-shape algorithm was proposed to construct a Delaunay triangulation network and obtain marine lane boundary. The theoretical analysis and experimental results show that the proposed method can achieve an accuracy of 86.7% and a recall rate of 79.4% when the maximum density value is 200， minimum density value is 10， length and width of window are 5 and 5 respectively. The experimental results show that the proposed method is effective to extract valuable marine lane boundaries from large-scale trajectory data.

Key words： trajectory data; Automatic Identification System （AIS）; spatio-temporal big data; Delaunay triangulation network; marine lane extraction

0 引言

由于傳统道路数据依靠人工测量、高分遥感影像等采集方式，具有价格昂贵、更新慢等特点，使得道路信息的快速获取和更新成为一直以来亟待解决和广受关注的问题。传统道路信息提取的难点在于人工测量及高分遥感影像的获取昂贵，以及基于图像的准确道路检测与识别技术。其中，人工测量在采集信息时耗时耗力，而且受到地理环境和自然条件等限制，会带来安全问题和时间成本大问题。高分遥感影像通过卫星进行高清影像数据采集，需要对数以百万量不同级别的移动物体长期且频繁地进行数据采集和处理，由于成本高昂难以推广使用。尤其是，基于图像的识别检测技术对已建设的陆地道路路网识别是一种常用的方式，但海洋航道不同于陆地道路，该方法无法基于海洋影像识别航道与非航道区域。近年来，随着移动传感器、物联网、云计算、大数据等技术的发展，产生了海量的时空轨迹大数据，并能够对其进行处理和分析。通过对众多在道路上行驶的移动物体产生的时空轨迹大数据进行分析，提取其中的道路信息成为一种受关注的技术。与陆地上的道路相比，海上没有人为建筑的航道，且受气候、洋流、事故等各种外部因素影响动态变化，情形更加复杂。与陆地道路信息采集相比，利用传统的人工测量等方法采集航道信息代价更为高昂，而利用来自船舶自动识别系统（Automatic Identification System， AIS）的众源船舶轨迹数据提取航道的信息，如果可行其一定具有低廉、更新快等特性，不失为一种可以探索且较有前景的途径;但是，众源船舶轨迹数据具有海量、地理范围大、数据质量差、数据密度在不同区域差异明显等特征[1]，这给海上航道的获取和更新提出了更大的挑战，本文提出了基于众源船舶轨迹数据进行航道信息提取的一整套方法，所瞄准的主要挑战及本文贡献包括：

1）针对众源船舶轨迹数据规模大的特点，采用基于MapReduce的方法对众源船舶轨迹数据进行并行化预处理，并基于并行化的Geohash编码方法进行空间聚类，从而既保留了船舶轨迹数据中隐藏的航道信息，又将海量的众源船舶轨迹数据进行了有效化简[2-3]。

2）针对大范围众源船舶轨迹数据在不同区域分布不均的问题，提出SpatialNiBlack局部动态阈值过滤算法[4]，对轨迹数据的聚类结果进行局部阈值过滤，并提出航道边界提取算法del-alpha-shape，能够对大范围的众源轨迹数据进行航道平面边界的准确提取[5]。

需要指出的是，航道本是具有寬度、水深等属性的水域，本文主要关注航道的平面属性，聚焦于如何基于海量的船舶轨迹数据提取航道的平面边界信息。

1 相关工作

当前，利用船舶轨迹数据进行航道提取的同类研究还比较前沿，直接相关的工作还比较少。相关工作有两类：一类工作通过对轨迹进行行为模式的统计，基于移动对象的行为模式与位置对比分析进一步建立航道[6-10];另外一类是利用K-Means、DBSCAN（Density-Based Spatial Clustering of Applications with Noise）等聚类算法挖掘航道的行为模式进一步提取航道[11-12]。虽然利用船舶轨迹数据进行航道提取的工作较少，但利用众源轨迹进行道路边界的提取是交通大数据领域的一个逐渐受到越来越多的研究者关注的研究问题，国内外学者主要围绕道路边界的提取和功能区域边界的提取展开了一定的研究工作，已经有一些值得借鉴的相关工作[13-14]。这些工作主要通过约束构造平面点集形状实现边界提取。如文献[15]运用约束德洛内（Delaunay）三角网从车辆轨迹线集中提取道路边界[15]，主要通过三角形边长度和Voronoi面积等几何特征表达轨迹点分布的聚集性差异，然后将两种不同维数的控制条件集成建立道路边界识别模型，运用“种子点”扩展方法实现道路边界的提取。也有一些学者提出双阈值的Alpha Shapes算法[16]，利用简单环的概念设计轮廓搜索算法，获得既有较好完整性又有较高几何精度的建筑物轮廓线。有的学者是基于凸壳的方法，采用凸壳的内缩特征获取边界的凹凸特征，这种方法一定程度上能处理空间密度分布不均的情况，但难以保证边界的规则性，也无法识别空洞现象[17]。从内容上看，这些方法无法适应采样稀疏、高噪声的全球定位系统（Global Positioning System， GPS）船舶轨迹数据，针对大规模船舶轨迹数据难以取得良好的效果。

综上分析，现有的平面点集形状重构算法能够通过多条件约束表征点集的形状，但是对于大规模的船舶轨迹数据则存在以下两点问题：1）由于采样环境原因，原始的大规模船舶轨迹数据包含大量噪声数据和丢失数据，增加了边界提取的难度;2）海洋上的轨迹数据受船舶的行驶习惯和活动热度影响，不同区域的船只由于活动程度不同，造成同一时间范围下轨迹数据量的量级在不同区域存在差异，使提取航道边界问题变得更加复杂。针对以上存在的问题，本文提出了一种新的算法模型，首先针对大规模的船舶轨迹数据，使用并行技术进行插值、去噪和Geohash聚类预处理;然后通过时空NiBlack局部动态阈值过滤和约束三角构网实现航道边界的提取。

2 方法概述

本文提出的算法是一种运用空间NiBlack局部动态阈值过滤及约束三角网从海量船舶轨迹点集中提取航道边界的过程。首先，对海量的船舶轨迹数据进行并行的去除异常轨迹点以及缺失数据补全预处理，提升数据质量;然后利用并行的Geohash编码方法对其进行聚类，以在保障结果精度的前提下提高后续数据处理和分析的效率;其次，提出空间NiBlack算法，针对大范围海量船舶轨迹数据在不同区域分布不均的问题，对聚类结果进行局部阈值过滤，对轨迹数据的预处理结果进行航道识别;最后，基于alpha shape边界算法提出一种新的提取算法del-alpha-shape，利用上一步局部阈值过滤的结果基于Delaunay三角剖分算法构造三角网，并提取边界轮廓，经过阈值过滤等处理后得到航道平面边界。

3 并行化轨迹数据空间聚类

3.1 并行化预处理

由于测量误差、传感器误差等因素，GPS轨迹数据存在数据缺失和噪声的问题，这样的数据在航道提取中很难发挥作用，所以在应用中需要对原始数据作预处理。本文轨迹预处理过程主要由过滤和插值两部分组成，过滤是为了消除噪声，插值是为了补全缺失数据。因为船舶在正常航行过程中不会频繁转向，产生的轨迹在一定范围内近似于直线，所以为了提高程序效率，本文采用基于时间和空间阈值的过滤和插值算法，而不是其他复杂算法。过滤算法的原理是计算点到轨迹的距离，如果点到轨迹的距离大于预定阈值则删除该点，从而消除孤立的噪声点。插值算法的原理是计算两个相邻轨迹点的距离和整体轨迹相邻两点的平均距离，根据相邻两点的距离与平均距离的比值大小，在这两个点之间插入不同数量的轨迹点，比值越大插入的点就越多。

本文采用某海洋船舶监测系统所采集的12个月的数据作为数据源，其中一个月约有60GB的数据，3亿多条记录。在使用CentOS 6.7单机实验时，对于60GB的数据量预处理完成大概需要40h，因此需要采用基于集群的并行化处理技术进行效率的提升。这些数据以SequenceFile序列化文件类型保存在Hadoop分布式文件系统（Hadoop Distributed File System， HDFS）中。因为原始数据集中每条记录均包含多个字段，本文对原始数据进行筛选，每条记录处理为定义1的四个基本字段组成。

本文方法基于MapReduce对数据并行化预处理，算法流程如图1所示。

1）数据分割，将HDFS中保存原始数据的SequenceFile文件划分成m（mn，n表示数据节点数）个数据块Split，每个块交由一个数据节点Datanode处理。

2）Map阶段，程序逐行读取字段无缺失数据，提取每条数据的v、x、y、t四个属性，以字段v为键，元组（x，y，t）为键值输出，其形式为〈v，（x，y，t）〉。

3）在Reduce阶段，每个Reduce处理具有相同键v的数据，操作如下：

第1步 将相同v的数据按时间t排序，设定时间阈值t_z和距离阈值d_z，轨迹阈值n。

第2步 将第i个轨迹点保存在数组list中，如果第i+1个轨迹点和第i个轨迹点之间的距离di小于d_z，且时间间隔Δti

第3步 計算数组list中轨迹点的个数N，如果N

第4步 计算相邻两点第j和第j+1之间的距离dj，如果dj小于d，或者两点的经度之差大于300（分别在东经180°附近和西经180°附近的相邻两点之间不用插值，这个判断的值可以是小于360°且尽可能大的值，设为300），则j=j+1，继续执行此步;否则令s=dj/d+2，在第j和第j+1两点之间插入s个点，并保存在list中， j=j+1，继续执行此步，直到遍历list，执行第5步。

第5步 保存list中的元素，并清空list，i=i+1，执行第2步，直到遍历相同v的数据。将预处理后的数据以字段v为键、元组（x，y，t）为键值输出，其形式为〈v，（x，y，z）〉。

4）将预处理后的数据按v保存在不同的文本文件中，最后得到预处理后的数据。在使用上述基于MapReduce的数据预处理之后，在表1所述的集群实验条件下，整个过程从40h缩短为150min。

3.2 基于Geohash编码的并行化空间聚类算法

轨迹聚类采用基于Geohash编码的算法，当轨迹点较多时，对于轨迹点，并不需要遍历计算数据集中所有轨迹点与当前点的距离，仅需关注小范围内的相邻轨迹点，因此考虑引入Geohash空间编码技术并行化聚类GPS轨迹点。基于Geohash编码空间聚类算法的原理是将地球看作一个平面，把地球平面划分为大小相同的网格，所有的轨迹点都会对应一个栅格，从而将位于同一栅格的轨迹点聚类。在使用CentOS 6.7单机实验时，对于60GB的数据量Geohash编码完成大概需要18h，因此需要基于集群进行并行化处理。

本方法基于MapReduce对数据并行处理，算法流程如图2所示。

程序的输入数据是3.1节预处理后的数据，程序按照以下流程运行。

1）数据分割，将HDFS中以不同v命名的m个数据文件划分成m′（m′≥n，n表示数据节点数）个数据块Split，每个块交由一个Datanode处理;

2）Map阶段，抽取每条数据的经度x和纬度y并进行Geohash编码，然后以编码code为键，1为键值输出，其形式为〈code，1〉;

3）Reduce阶段，统计不同Geohash编码code的密度dsy，计算其对应区域C的中心点经纬度C_x和C_y，以C_x、C_y和dsy为键，null为键值输出，其形式为〈（C_x，C_y，dsy），null〉;

4）将所有数据按元组〈C_x，C_y，dsy〉保存。完成空间聚类。

在本文中，Geohash编码位数为11位，大约等于11km见方的一个区域。在使用上述基于MapReduce的数据预处理之后，在表1所述的集群实验条件下，整个过程从18h缩短为110min。

4 基于NiBlack局部动态阈值过滤及alpha shape的边界提取算法

4.1 SpatialNiBlack局部动态阈值过滤算法

Niblack算法是一种比较常用、简单、有效的局部动态阈值算法，通常应用于数字图像分割、图像二值化等领域，该算法根据局部平均值和局部标准差确定当前区域阈值。由于不同区域的轨迹点密度存在差异，固定阈值无法准确地过滤数据，因此本文基于NiBlack动态窗口过滤算法思想，提出适于轨迹数据分析的SpatialNiBlack算法，实现不同密度的轨迹点数据过滤，该方法有效克服了固定阈值的缺陷。

本文所用的数据是3.2节中聚类的结果，这些数据被保存在一个二维数组中，数组中的每一行用一个三元组〈C_x，C_y，dsy〉来表示一个栅格C，其中C_x和C_y分别表示栅格中心点的经度和纬度，dsy表示栅格C中轨迹点的密度。

本方法步骤如图3所示，基于SpatialNiBlack局部动态阈值的过滤算法如下。

1）首先得到全球11位的栅格，设置大小为m×n的矩形窗口，因为数据从纬度变化方向来看，有密度差异，所以窗口设置成扁平的矩形，使窗口内大多的点的密度在同一个级别。窗口包含m×n个11位栅格，执行第2）步;

2）遍历当前窗口中的栅格，统计窗口中所有栅格的dsy值，根据式（1）计算出阈值T（x，y），然后判断窗口中心点所在栅格是否保留。执行下一步;

其中对于栅格中心点经纬度（x，y），T（x，y）为该位置的阈值，m（x，y）为该点所在窗口内轨迹点密度的均值，s（x，y）为该点所在窗口内轨迹点密度的标准差，k为修正系数（通常取-0.1）。

3）窗口中心栅格后移一个格，以该栅格为中心建立新的窗口，执行第2）步;直到窗口遍历所有栅格。最后输出过滤后的数据。

4.2 基于alpha shape的边界提取

本文基于alpha shape边界算法，提出一种新的提取算法del-alpha-shape，该算法能够有效地对点集进行边界提取。

alpha shape是一种经典的点集轮廓提取算法，即输入点集，通过对任意两点建立半径为alpha的外接圆，如果没有第三个点在该圆内，则这两点连成的线段是边界线段。通过对alpha值的调整可以得到不同精确程度的轮廓。如果两个点连线距离大于2倍alpha阈值，那么无法构成半径alpha的外接圆，对这两点不进行判断。

Delaunay三角剖分算法是将点集处理成由三角面构成的凸包平面图，该平面图内除了端点，不包含任何点集中的点。

定义9 三角形密度指数alpha。三角形的三个点可以构成唯一的一个外接圆，该外接圆半径circum_r的倒数值，稱为三角形密度指数。

在Delaunay三角网图（图4）中，航道内的三角形的密度指数较大，航道外的三角形的密度指数较小。设定alpha阈值（用alpha_value表示），如果三角形alpha≥alpha_value，为航道三角形，保留该三角形的三条边到Edges;否则，为空洞三角形，不保留边。

该算法流程如图5所示。

本文4.1节得到的栅格过滤结果是一个点集数据，以Points表示。该算法首先使用Delaunay方法对点集Points进行Delaunay三角化，得到三角面集合Triangles（第1）行）;然后遍历Triangles，使用calCircum方法计算每个三角形的密度指数circum，如果circum不小于密度指数阈值，则将该三角形的边加入边集Edges（第2）～6）行）;下一步用polygonize方法对Edges进行多边形化，构成多边形集Polys（第7）行）;最后，遍历Polys，判断每个多边形的面积是否大于面积阈值maxArea，如果大于，则将该多边形的边界坐标polygon.coords添加到Channels集合（第8）～11）行）。算法时间成本主要在第7）行多边形生成过程，整体时间复杂度是O（N2）。

5 实验分析

5.1 数据集介绍

船只航行轨迹数据是全球2016年内的所有船只轨迹。船只航行轨迹数据包括GPS时间戳、水上移动通信业务标识码（MMSI）、GPS经度、GPS纬度等多项信息。轨迹点采样间隔在10～120s，总数据量在1TB左右。

本实验选取两次不同的数据：一是2016年6月渤海区域的油轮航行轨迹数据;一是2016年6、7、8月“一带一路”区域大宗货物船只航行轨迹数据。本实验的预处理步骤是在6台搭建了Hadoop环境的机器下进行的，系统为CentOS release 7.0操作系统，JDK版本为1.7.0，Hadoop版本为2.6.0。集群具体配置如表1所示。

实验合并过滤程序使用Java1.8编程语言实现，航道提取使用Python 3.5编程语言实现，在单台物理机（64GB内存，系统CentOS release 7.0）运行。

5.2 实验结果与分析

本文设计以下实验验证模型，第一组实验使用数据是中国—东南亚—非洲连线区域2016年2—12月的大宗船只的航行数据。

图6（a）是使用3.1和3.2节介绍预处理方法，预处理相关的参数取值为：距离阈值30km，时间阈值3600s，分段轨迹点数目阈值30，对中国—东南亚—非洲连线区域2016年2—12月的大宗船只原始轨迹数据预处理后的栅格化聚类结果上图此句不通顺，请作相应调整。。图6（b）是使用4.1节介绍方法滑动窗口局部动态阈值过滤的结果。这里涉及到几个参数设置，栅格位数11位、栅格最大阈值200、窗口宽50、窗口高3。图6（c）是基于滑动窗口局部动态阈值过滤结果用4.2节介绍方法作航道提取的上图图6（c）是基于滑动窗口局部动态阈值过滤和4.2节介绍方法做航道提取的结果[18]此句不通顺，请作相应调整。。

第二组实验使用数据是渤海区域的2016年6—8月的油轮航行数据，主要为了进一步验证较小地理空间范围下模型提取航道和空洞的准确性。

图7（a）是使用3.1和3.2节介绍的预处理方法，预处理相关的参数取值为：距离阈值30km，时间阈值3600s，分段轨迹点数目阈值30，对渤海区域2016年6—8月油轮原始轨迹数据预处理后的栅格化聚类结果如图7。

图7（b）是使用4.1节介绍方法滑动窗口局部动态阈值过滤的结果。这里参数设置为，栅格位数16位，栅格最大阈值500，窗口宽15，窗口高7。图7（c）是基于滑动窗口局部动态阈值过滤结果用4.2节介绍方法作航道提取的上图图7（c）是基于滑动窗口局部动态阈值过滤和4.2节介绍方法做航道提取的结果。

本文以渤海某一区域的使用固定阈值进行航道提取的结果为参考，如图8（a），然后使用本文提出的大规模窗口局部过滤法的航道结果与之对比分析。采用过滤结果叠置方法评价航道提取的有效性，并借鉴文献的方法计算两种常用的评价指标，即准确率（precision）和召回率（recall）来评价叠置结果。准确率反映实验结果的准确度，召回率反映实验结果的完整度。

本文准确率和召回率的计算如式（2）和式（3）：

其中：gridNumextstd表示新提取的航道栅格结果与标准航道栅格结果重合的数目，gridNumstd表示标准航道栅格结果的数目，gridNumext表示提取航道栅格结果的数目。

实验数据是2016年6月至2016年—8月份油轮轨迹数据，在渤海某一区域内进行航道提取方法验证，采取不同参数对比结果如图8。

参数值含义

maxDsy100～∞这个数据项是何意？100后面加了“--”，何意？回复：表示取值范围，可改为100-∞最大密度值

minDsy0～10最小密度值

w_size3×3、5×5、10×3滑动窗口大小

表3是选用不同参数进行多组实验后，与标准航道结果叠置后求得的准确率和召回率结果。其中：precision表示准确率，recall表示召回率。从表中可以看出：

1）最大密度值maxDsy的限制能够避免局部窗口内方差过大造成的阈值过大，图8（b）导致窗口内所有值被过滤舍弃，造成召回率偏低，因此，设置最大密度值进行限制，并分多组实验测试当该值为200时效果准确率和召回率都有较好的表现。

2）最低密度值minDsy可以过滤掉一些存在的噪声数据，提升召回率，经过多组测试发现10是一个比较合适的值。

3）滑动窗口的大小选择对结果有一定影响，当窗口较小时，准确率因阈值计算考虑范围过小造成准确率下降，如图8（k）;当窗口形状不合适，可能造成阈值结果过大，过滤舍弃掉过多的值导致航道部分缺失，如图8（l）。

实验验证了本文所提方法能够在大范围数据下保持较高的提取精确度，同时依然保持较高的完整度;同时，本文所提算法与其他相似的研究算法有所区别：

1）基于三角网的边界提取。文献[15]通过对轨迹数据构造三角网并约束条件提取道路，但其主要思路是通过自定义的长短边概念找到边界。本文主要通过三角网的外接圆半径约束条件寻找边界边;同时，通过实现两种算法进行实验，发现本文算法在相同量的大规模数据下的效率更高，运行时间效率可提高约5倍此处表述不当，意义不对，“运行时间”应该是“减少几分之几”，请按照这个思路进行调整;或者改为“运行效率可提高约5倍”，这样修改合适吗？请明确。

2）基于GPS道路轨迹线提取。文献[13]通过聚类方式改进算法，提取了不同级别道路的中心线结果。本文所提航道除了能够表示路网结构外，还能够看出航道的边界大致范围，由于航道在海面上没有宽度限制，因此该范围是一个不严格的大概范围。

3）图像轮廓线提取算法。文献[16]借助图像处理的研究中进行，其输入大致是图像数据，通过对像素值的识别计算，提取出包括空洞的轮廓线结果。本文算法是基于时空轨迹大数据进行的边界提取，由于时空轨迹数据密度值相比0～255的像素值范围具有更大的取值范围，因此在算法思路上有所区别，本文在4.1节局部过滤算法内通过参数的限制有效解决了时空数据复杂的问题。

6 结语

本文提出一种基于大规模船舶轨迹数据进行航道边界提取的方法，通过去噪、插值、分段预处理和Geohash编码简化处理数据;然后基于处理后的数据结果使用文中提出的SpatialNiBlack航道识别算法过滤航道;最后使用del-alpha-shape算法获得有效的航道边界。通过实验可以发现：该方法能够有效解决边界提取的空洞问题;在较大的地理范围下能够进行航道边界提取;通过改进的局部过滤算法，可以解决全局范围下的数据因存在质量差异导致部分航道丢失的问题。进一步改进的方面：完善算法，使其能夠同时对近海区域的一些较为精细的航道进行更加有效的提取。

参考文献 （References）

[1] WANG Y， ZHU Y， HE Z， et al. Challenges and opportunities in exploiting large-scale GPS probe data， HPL-2011-109 [R]. Palo Alto， CA： HP Laboratories， 2011.

[2] ZHANG W， GOERLANDT F， KUJALA P， et al. An advanced method for detecting possible near miss ship collisions from AIS data[J]. Ocean Engineering， 2015， 107（1）：60-69.

[3] WANG X， LIU X， LIU B， et al. Vessel route anomaly detection with Hadoop MapReduce[C]// Proceedings of the 2015 IEEE International Conference on Big Data. Piscataway， NJ： IEEE， 2015：25-30.

[4] NIBLACK W. An Introduction to Digital Image Processing[M]. Englewood Cliffs， NJ： Prentice-Hall， 1986：115-116.

[5] EDELSBRUNNER H， KIRKPATRICK D， SEIDEL R. On the shape of a set of points in the plane[J]. IEEE Transactions on Information Theory，1983，29（4） ：551-559.

[6] SUN F， DENG Y， DENG F， et al. Unsupervised maritime traffic pattern extraction from spatio-temporal data[C]// Proceedings of the 2015 International Conference on Natural Computation. Piscataway， NJ： IEEE， 2015：1218-1223.

[7] MAZZARELLA F， FERNANDEZ ARGUEDAS V， VESPE M. Knowledge-based vessel position prediction using historical AIS data[C]// Proceedings of the 2015 Sensor Data Fusion： Trends， Solutions， Applications. Piscataway， NJ： IEEE， 2015：1-6.

[8] FERNANDEZ ARGUEDAS V， PALLOTTA G， VESPE M. Maritime traffic networks： from historical positioning data to unsupervised maritime traffic monitoring[J]. IEEE Transactions on Intelligent Transportation Systems， 2017， 19（3）： 722-732.

[9] PALLOTTA G， VESPE M， BRYAN K. Vessel pattern knowledge discovery from AIS data： a framework for anomaly detection and route prediction[J]. Entropy， 2013， 15（6）：2218-2245.

[10] SPILIOPOULOS G， ZISSIS D， CHATZIKOKOLAKIS K. A big data driven approach to extracting global trade patterns[C]// Proceedings of the 2017 International Workshop on Mobility Analytics for Spatio-temporal and Social Data. Berlin： Springer， 2017：109-121.

[11] ETIENNE L， DEVOGELE T， BOUJU A. Spatio-temporal trajectory analysis of mobile objects following the same itinerary[C]// Proceedings of the 2010 International Symposium on Spatial Data Handling. Hong Kong： [s.n.]， 2010：86-91.

[12] 杨杰，李小平，陈湉.基于增量时空轨迹大数据的群体挖掘方法[J].计算机研究与发展，2014，51（s2）：76-85.（YANG J， LI X P， CHEN T. A group mining method for incremental spatio-temporal trajectory big data[J]. Journal of Computer Research and Development， 2014， 51（s2）：76-85.）

[13] 欧阳鸿，刘建勋，刘毅志，等.基于步行GPS轨迹的路网提取方法[J].計算机与现代化，2014（2）：124-128.（OUYANG H， LIU J X， LIU Y Z， et al. An extraction method of road network based on walking GPS trajectories[J]. Computer and Modernization， 2014，（2）：124-128.）

[14] 杨伟，艾廷华.基于车辆轨迹大数据的道路网更新方法研究[J].计算机研究与发展，2016，53（12）：2681-2693.（YANG W， AI T H. A method for road network updating based on vehicle trajectory big data[J]. Journal of Computer Research and Development， 2016，53（12）：2681-2693.）

[15] 杨伟，艾廷华.运用约束Delaunay三角网从众源轨迹线提取道路边界[J].测绘学报，2017，46（2）：237-245.（YANG W， AI T H. The extraction of road boundary from crowdsourcing trajectory using constrained Delaunay triangulation[J]. Acta Geodaetica et Cartographica Sinica， 2017， 46（2）： 237-245.）

[16] 李云帆，谭德宝，高广，等.双阈值Alpha Shapes算法提取点云建筑物轮廓研究[J].长江科学院院报，2016，33（11）：1-4.（LI Y F， TAN D B， GAO G， et al. Extraction of building contour from point clouds using dual threshold Alpha Shapes algorithm[J]. Journal of Yangtze River Scientific Research Institute， 2016，33（11）： 1-4.）

[17] 朱杰，孙毅中.多约束的平面点集形状重构方法[J].测绘学报，2017，46（2）：253-264.（ZHU J， SUN Y Z. An efficient approach to shape reconstruction from planar point set based on multi-constraints[J]. Acta Geodaetica et Cartographica Sinica， 2017，46（2）：253-264.）

[18] WU L， XU Y， WANG Q， et al. Mapping global shipping density from AIS data[J]. Journal of Navigation， 2017， 70（1）：67-81.