CN111291776B - 基于众源轨迹数据的航道信息提取方法 - Google Patents

Abstract

本发明提供一种基于众源轨迹数据的航道信息提取方法。该方法包括将指定的地理范围按照要求的精度划分为多个栅格，众源轨迹数据的每个轨迹点对应到相应的一个栅格，将位于同一个栅格的轨迹点使用该栅格的指定点来表示，获得化简后的数据；根据所述化简后的数据的轨迹点在栅格内的密集程度，按照设定的阈值对栅格进行合并，获得合并后的栅格集；从所述合并后的栅格集提取航道区域边界信息。利用本发明的方法能够快速准确的提取航道边界。

Description

基于众源轨迹数据的航道信息提取方法

技术领域

本发明涉及数据挖掘技术领域，尤其涉及一种基于众源轨迹数据的航道信息提取方法。

背景技术

道路相关地理数据是国家基础地理信息、智能交通的重要组成部分，在智慧城市建设、智能导航、交通管制、网络地图服务等方面具有重要应用价值。随着移动传感器、云计算等技术的发展，近年来，在交通和GIS相关应用领域，利用来自于大量交通工具(例如汽车、船舶等)的海量轨迹数据(也称为众源轨迹数据)提取道路的地理信息，相较于传统的道路地理信息获取方式，具有低廉、更新快等特性。通常，在城市交通领域，可以进行道路地理信息提取的众源轨迹数据采集自陆地交通工具的全球定位系统(GPS)终端设备或GPS采集中心，在海上交通领域，众源轨迹数据采集自船舶自动识别系统(AIS)终端设备或岸基AIS数据采集中心。

相对于陆地交通工具的GPS数据，基于船舶AIS数据的航道提取面临更大的挑战。这是因为：船舶众源轨迹数据具有规模大、噪声高、采样频率分布不均、以及密度、质量分布不均等特点，例如，全球船舶轨迹数据一年采集的原始数据为T级别，几乎每一只船舶每次出行的轨迹都存在错误的采样，近海区域轨迹点采样频率为5秒到100秒不等，远海区域轨迹点采样频率2分钟到10分钟不等，采样间隔较大；此外，船舶AIS数据来源各异，密度、质量分布不均。

在现有技术中，基于众源轨迹数据进行道路(或航道)信息提取的方法主要包括：1)、聚焦于用轨迹数据聚类的方法提取道路中心线，例如将K-Means聚类与高斯模型相结合提取道路中心线结构并识别车道；通过合并聚类轨迹线提取非交叉口道路中心线。但是聚类不适用大范围下、轨迹点稀疏的数据。2)、将航道看作多边形几何形状，利用众源轨迹数据的几何特征基于三角网等方法提取航道的边界，例如，运用约束Delaunay三角网从车辆轨迹线集中提取道路边界，或者将轨迹栅格化，利用矢量化算法提取道路面。但是现有基于三角网的方法无法处理海量的船舶轨迹数据。3)、利用图像处理技术提取道路骨架线并构建地图，例如，将轨迹点转为二值图像，利用形态学方法提取道路骨架线构架地图，或者将核密度估计与聚类结合提取路网信息。

目前的道路信息提取主要的问题有：1)、大多算法仅提取了道路结构中心线，没有准确地提取道路内外部边界信息。2)、大多研究针对陆地一定范围的轨迹数据进行边界提取，不适用于较大范围采样数据密度不均匀的海量数据。对海洋船舶众源轨迹数据而言，航道在远海和近海范围内采集得到的数据密度、数据质量差异很大，且近海区域的船舶轨迹点天然比远海区域的船舶轨迹点分布更为分散，对提取的精细程度要求也更高，但在远海区域船舶轨迹点分布较为集中，对航道提取的精细程度要求较低。总之，现有方法无法使用统一的精度范围去表示提取的航道，不适用于基于大范围船舶众源轨迹数据的航道提取。

因此，需要对现有技术进行改进，以提供更精确的航道提取方法。

发明内容

本发明的目的在于克服上述现有技术的缺陷，提供一种基于众源轨迹数据的航道信息提取方法。

根据本发明的第一方面，提供一种基于众源轨迹数据的航道信息提取方法。该方法包括：

步骤1：将指定的地理范围按照要求的精度划分为多个栅格，众源轨迹数据的每个轨迹点对应到相应的一个栅格，将位于同一个栅格的轨迹点使用该栅格的指定点来表示，获得化简后的数据；

步骤2：根据所述化简后的数据的轨迹点在栅格内的密集程度，按照设定的阈值对栅格进行合并，获得合并后的栅格集；

步骤3：从所述合并后的栅格集提取航道区域边界信息。

在一个实施例中，步骤1包括：对于所述多个栅格进行GeoHash编码，统计不同GeoHash编码的轨迹点密度，计算其对应栅格的中心点的经度和维度，并保存每个栅格的中心点的经度、维度和轨迹点密度作为化简后的数据。

在一个实施例中，步骤2包括：

步骤21：将获得的栅格数据利用四叉树结构存储，除根节点之外，每个节点存储一个栅格的GeoHash编码和轨迹点密度；

步骤22：对于四叉树结构中的父节点对应的四个子节点，在该四个子节点的轨迹点密度均小于设定的第一密度阈值时，合并该四个子节点。

在一个实施例中，在步骤22之后还包括：

步骤23：设定包含固定栅格个数的矩形作为窗口，以窗口边沿着整个地理范围边对齐后的窗口中心点形成的范围作为窗口中心点的滑动范围进行遍历，如果窗口中所有栅格的轨迹点密度大于设定的第二密度阈值，则保留该窗口内的栅格。

在一个实施例中，所述第二密度阈值设置为：

T＝avg+alpha1×var

其中，avg表示窗口内所有栅格的轨迹点密度的平均值，var表示窗口内所有栅格的轨迹点密度的方差，alpha1表示方差修正系数。

在一个实施例中，步骤3包括：

步骤31：利用Delaunay对合并后的数据进行三角化，得到三角面集合；

步骤32：计算三角面集合中每个三角形的密度指数，将密度指数小于指数阈值的三角形的边加入到边集；

步骤33：对边集进行多边形化，获得多边形集；

步骤34：将多边形集中面积大于预定的面积阈值的多边形的边界坐标作为航道边界的顶点坐标集合，以形成多边形的航道边界。

在一个实施例中，所述三角形的密度指数为：

alpha2＝1/circum_r

其中，circum_r表示三角形的外接圆半径。

在一个实施例中，在步骤1中，所述众源轨迹数据是经过预处理之后的获得的数据，预处理过程包括：

排序步骤：以时间为顺序对众源轨迹数据进行排序；和/或

采样步骤：基于众源轨迹数据的相邻两个轨迹点的时间间隔与第一时间间隔阈值的比较结果确定需要保留的轨迹点；和/或

过滤步骤：根据众源轨迹数据的每个轨迹点的速度与速度阈值的比较结果确定需要过滤掉的轨迹点；和/或

插值步骤：根据过滤后的相邻轨迹点之间的时间间隔与第二时间间隔阈值的比较结果确定在相邻轨迹点之间需要插入的轨迹点。

在一个实施例中，通过以下步骤执行所述预处理过程：

数据分割步骤：将原始众源轨迹数据划分成m个数据块，每个数据块交由一个数据节点处理，其中，m＞＞n，n表示数据节点数；

Map阶段：对于所分配的数据块，多个数据节点并行逐行读取字段无缺失数据，提取每条数据的v、x、y、t四个属性，以字段v为键，元组(x，y，t)为键值输出，其中，v是船舶的唯一标识，x是轨迹点的经度，y是轨迹点的维度，t是时间戳；

Reduce阶段：每个Reduce处理具有相同键v的数据，执行所述排序步骤、采样步骤、过滤步骤和插值步骤。

与现有技术相比，本发明的优点在于：利用大规模的众源轨迹数据，针对包括远海和近海不同精细程度的大范围航道，通过并行化的自适应精度的合并过滤算法，能够进行不同精细程度的统一的航道识别及提取。

附图说明

以下附图仅对本发明作示意性的说明和解释，并不用于限定本发明的范围，其中：

图1示出了根据本发明一个实施例的基于众源轨迹数据的航道信息提取方法的流程图；

图2示出了根据本发明的一个实施例对众源轨迹数据进行采样、过滤和插值处理的示意图；

图3示出了根据本发明一个实施例的基于MapReduce的进行预处理过程的示意图；

图4示出了预处理过程的数据示例；

图5示出了预处理前后的轨迹效果对比图；

图6示出了根据本发明一个实施例的栅格化的过程示意图；

图7示出了根据本发明的一个实施例基于MapReduce对轨迹点进行化简的示意图；

图8示出了对东海区域栅格化的效果示意图；

图9示出了根据本发明的一个实施例建立四叉树的示意图；

图10示出了根据本发明一个实施例的基于四叉树对栅格合并的示意图；

图11示出了根据本发明一个实施例的栅格局部过滤的示意图；

图12示出了滑动窗口示意图；

图13示出了相邻栅格越界的示意图；

图14示出了滑动窗口中心范围的示意图；

图15示出了窗口中心点遍历的示意图；

图16示出了东海区域合并过滤后的效果示意图；

图17示出了根据本发明一个实施例的航道边界提取的流程图；

图18示出了某海域船舶轨迹的Delaunary的三角化图；

图19示出了三角化示意图；

图20示出了多角化示意图；

图21示出了多角化过滤的效果图；

图22示出了根据本发明的一个实施例对东海区域航道提取结果的示意图。

具体实施方式

为了使本发明的目的、技术方案、设计方法及优点更加清楚明了，以下结合附图通过具体实施例对本发明进一步详细说明。应当理解，此处所描述的具体实施例仅用于解释本发明，并不用于限定本发明。

为了便于理解，首先介绍本发明涉及的一些概念：航道，将航道定义为一个二维多边形形状，代表海上船舶在规范下能够航行的区域范围；空洞，是由于航道内部存在礁石等障碍物，而出现的不能被船舶行驶的区域范围；非航道，除航道和空洞以外的船舶不允许行驶的区域范围。栅格，是指地图上的一个矩形区域，通过对地理空间进行网格划分，把空间范围划分为大小相同的网格，每一个网格称为栅格，栅格编码的位数称为栅格精度，栅格精度越高，对应的栅格面积越小。而一个栅格内包含的轨迹点的数目称为栅格密度；航道边界，航道边界可精确表示，航道平面边界是平面多边形的集合，表示为C＝{c₁,c₂,…,c_i,…,c_n}，其中，c_i＝{<x_i1,y_i1>,<x_i2,y_i2>,…,<x_in,y_in>}是一个多边形，x_ij，y_ij分别表示第i个多边形的第j个顶点的经度坐标和纬度坐标,多边形c_i可由其顶点集合按照顺时针或逆时针方向的序列来表示；航道精度，提取航道所用栅格的栅格精度的均值称为航道精度。

图1示出了根据本发明一个实施例的基于众源轨迹数据的航道信息提取方法，简言之，该方法包括对众源轨迹数据进行预处理、利用栅格化和编码方式简化轨迹点、对轨迹点较稀疏的栅格进行合并、过滤非航道栅格和提取多边形航道边界等，从而获得航道边界信息，具体地，包括以下步骤：

步骤S110，对众源轨迹数据进行预处理。

预处理的过程可包括对众源轨迹数据进行排序、采样、去噪过滤和插值等中的一项或多项，以去除噪声和存在信息丢失的数据，而保留对提取航道信息有重要意义的数据。例如，排序可将每条船的轨迹点数据按时间先后顺序进行排列；采样可对时间间隔过于小的轨迹点进行过滤；去噪可对船舶速度不符合常理、包含轨迹点过少的轨迹段进行过滤；插值可补齐部分缺失数据。

下文将分别介绍预处理过程的实施例。

1)、排序实施例

通过排序可将众源轨迹数据按照时间顺序进行排列，以便于后续处理。

例如，采用某海洋船舶监测系统所采集的12个月的数据作为数据源，一个月约有60GB的数据，3亿多条记录。这些数据可以SequenceFile文件类型保存在Hadoop的分布式文件系统HDFS中。由于原始数据集中每条记录均包含多个字段，在对原始数据进行初步筛选后，保留每条记录的时间戳、船舶mmsi号(即船舶的唯一标识)、经度和纬度四个基本字段，对原始数据进行初步筛选后的样例数据形式如下：

**********

1456714928，100900074，119.1305，39.091317

1456754719，100900074，119.174，38.922283

************

基于时间戳进行排序之后，可获得按照时间顺序排列的数据。

2)、采样实施例

在一个实施例中，通过采样将对时间间隔小于阈值的轨迹点进行过滤。

参见图2(a)所示，原始轨迹L由序列{P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13,P14,P15}组成，下表1中记录了原始轨迹中每一个轨迹点的时刻和与上一个轨迹点之间的距离。

表1：轨迹点时刻和各轨迹点间的距离

假设采样的时间间隔阈值设置为10s，计算相邻两个轨迹点之间的时间间隔，如果不小于10s则保留，否则删除。具体地，从P1开始，计算P1和P2的时间间隔为5s，小于10s，因此删除P2点，然后再计算P1和P3的时间间隔为10s，不小于10s，因此保留P3点；接着从P3点开始，和后面的点比较，每次删除一个点后，则用这两个点中的前一个点和删除点的后一个点比较，每次保留一个点后，则从保留的点开始和后面的点比较。以此类推，直到遍历所有轨迹点。通过这种方式，经过采样后的轨迹L由序列{P1,P3,P5,P6,P7,P8,P9,P10,P11,P12,P14,P15}组成，参见图2(b)所示。

3)去噪过滤实施例

通过去噪可以对船舶速度不符合常理、包含轨迹点过少的轨迹段进行过滤。

在一个实施例中，用相邻两点之间的距离除以相邻两点的时间间隔作为后一个点的速度，计算轨迹中每个点的速度，如果速度大于速度阈值，则过滤掉该点。例如，以采样后的轨迹点为例，速度阈值设置为5m/s，计算采样后各轨迹点的速度，先保存P1点，用P1和P3点的距离和时间间隔计算P3点的速度，为1.8m/s，以此类推，P5点的速度为1.7m/s，P6点的速度为2m/s，P7点的速度为2m/s，P8点的速度为1.75m/s，P9点的速度为1.2m/s，P10点的速度为10m/s，由于P10点的速度大于设置的速度阈值，因此将点P10过滤掉，用P9和P11两点计算P11点的速度，为2m/s，然后计算P12点的速度为2m/s，P14点的速度为2.1m/s，P15点的速度为2.2m/s。通过这种方式，经过去噪过滤后的轨迹L由{P1,P3,P5,P6,P7,P8,P9,P11,P12,P14,P15}组成，参见图2(c)所示。

在此实施例中，采用了基于时间和空间阈值的去噪过滤方式，能够提高处理效率。

4)插值实施例

在一个实施例中，通过插值补齐缺失的数据，基本原理是计算两个相邻轨迹点的距离和整体轨迹上所有相邻两点的平均距离，根据相邻两点的距离与平均距离的比值大小，在这两个点之间插入不同数量的轨迹点，比值越大插入的点就越多。

在另一个实施例中，根据两个相邻轨迹点之间的时间间隔和采样阈值进行比较来确定是否在两个点之间插入轨迹点以及插入轨迹点的数目。

例如，计算过滤后相邻轨迹点之间的时间间隔，如果时间间隔t(i,i+1)大于采样阈值T_z(如10s)，则需要在这两个点之间线性插入n个轨迹点，其中n＝t(i,i+1)/T_z-1。具体地，从第一个点开始，计算P1和P3之间的时间间隔为10s，不需要插值；P3和P5之间的时间间隔为10s，不需要插值；P5和p6之间的时间间隔为30s，需要插入n个点，其中n＝30/10-1＝2；P6和P7之间的时间间隔为10s，不需要插值；P7和P8之间的时间间隔为20s，需插入n个点，其中n＝20/10-1＝1；P8和P9之间的时间间隔为10s，不需要插值；同理，遍历处理整个轨迹。通过这种方式，经过插值之后，轨迹L由序列由{P1,P3,P5,C1,C2,P6,P7,C3,P8,P9,P11,P12,P14,P15}组成，C1、C2、C3为插入的点，参见图2(d)所示。

本发明预处理过程的实施例，可采用并行分布式计算方法实现，以提高大规模数据量的处理效率。例如，采用基于MapReduce的算法，参见图3所示，其具体过程包括：

第一步、数据分割

将原始数据划分成m个数据块，每个块交由一个数据节点处理，其中，m＞＞k，k表示数据节点数。

第二步、Map阶段：并行提取数据

在Map阶段，逐行读取字段无缺失数据，提取每条数据的v、x、y、t四个属性，以字段v为键，元组(x，y，t)为键值输出，其形式为<v，(x，y，t)>，其中，v是船舶的唯一标识，x表示轨迹点的经度，y表示轨迹点的维度，t表示时间戳。

第三步、Reduce阶段：并行处理排序、采样、过滤和插值

在reduce阶段，每个Reduce处理具有相同键v的数据，先把数据按t排序，然后进行去噪过滤和插值等。例如，具体步骤如下：

步骤S111-a，将相同v的数据按时间t排序，设定采样间隔t_z、分段阈值t_s、速度阈值v_z和轨迹点阈值n；

分段阈值是指用于轨迹分段的相邻轨迹点时间阈值，当相邻的两个轨迹点之间的时间间隔小于该阈值时，可以将两个轨迹点看作属于同一次出行的两个轨迹点，即两个轨迹点属于同一个轨迹段。

步骤S111-b，对排好序的数据按采样间隔t_z采样，采样后相邻两个轨迹点的时间间隔不小于t_z；

步骤S111-c，将第i个轨迹点保存在数组list中，如果第i+1个轨迹点和第i个轨迹点之间的时间间隔Δt_i小于分段阈值t_s，则令i＝i+1(保存第i+1个点)，继续执行此步；否则，执行步骤S111-d；

步骤S111-d，计算数组list中轨迹点的个数N，如果N小于n(即一个轨迹段所包含的轨迹点太少)，则将数组list中的点视为噪点，并清空list，且i＝i+1，执行步骤S111-c；否则，将数组list中的点视为一段轨迹，执行步骤S111-e；

步骤S111-e，计算list中第j+1个点的速度v_j+1＝Dis(j,j+1)/T(j,j+1)，如果v_j+1小于速度阈值v_z，则保存第j+1个点；否则删除第j+1个点，直到遍历list；

步骤S111-f，计算相邻两点第j和第j+1之间的时间间隔Δt_j，如果Δt_j不大于t_z，或者两点的经度之差大于300(分别在东经180度附近和西经180度附近的相邻两点之间不用插值，这个判断的值可以是小于360度且尽可能大的值，设为300)，或者轨迹点的速度小于速度最小阈值，则令j＝j+1，继续执行此步；否则令s＝Δt_j/t_z，在第j和第j+1两点之间插入s-1个点，并保存在list中，令j＝j+1，继续执行此步，直到遍历list，执行步骤S111-g；

步骤S111-g，保存list中的元素，并清空list，i＝i+1，执行第三步，直到遍历相同v的数据。将预处理后的数据以字段v为键、元组(x，y，t)为键值输出，其形式为<v，(x，y，t)>。

图4示出了预处理过程的数据示例，其中，图4(a)是原始数据，图中数据的第一行是时间戳，可以看出，原始数据是无时间顺序的；图4(b)是排序后的数据，可以看出，相邻两点之间的时间间隔是不确定的，可能存在数据缺失的问题。图4(c)是经过采样、过滤和插值之后的数据，采样和插值的时间间隔是30秒，可见经过预处理之后，相邻两点之间的时间间隔是固定的，都是30s。

图5示出了预处理前后的轨迹效果对比图，图5(a)为预处理之前的轨迹，图5(b)是预处理之后的轨迹，可以看出，预处理之后的轨迹更加均匀和连续。

对于预处理之后的数据，可按照船舶的唯一标识，即按照v保存在不同的文本文件中。

步骤S120，通过栅格化对预处理后的数据进行化简，获得化简后的轨迹点。

当轨迹点较多时，不需要遍历计算数据集中所有轨迹点与当前点的距离，仅需关注小范围内的相邻轨迹点，因此可以对预处理之后的数据进行化简。

在一个实施例中，通过栅格化，采用Geohash空间编码技术对轨迹点进行化简，其原理是将地球看作一个平面，把地球平面划分为大小相同的网格(即栅格)，所有的轨迹点都会对应一个栅格，将位于同一栅格的轨迹点使用该栅格的中心点来代表。

图6示出了栅格化位数为1的栅格化过程，左图是地球平面，分布有一些轨迹点(黑色点)，首先，进行1次网格划分，将地球表面划分为4的1次方个网格P1、P2、P3、P4(如右图所示)，P1中心点坐标为(-90,45)(分别表示经度和维度)，P2的中心点坐标为(90,45)，P3的中心点坐标为(-90,-45)，P4的中心点坐标为(90,-45)；然后，统计每个区域轨迹点的个数(即每个栅格的密度)，P1的轨迹点密度dsy1为9，P2的轨迹点密度dsy2为3，P3的轨迹点密度dsy3为9，P4的轨迹点密度dsy4为5；最后，保存栅格的中心点坐标和轨迹点的密度，即，P1(-90，45，9)；P2(90，45，3)；P3(-90，-45，9)；P4(90，-45，5)。类似地，可获得其他精度的栅格化的数据。

由于众源轨迹数据规模较大，例如，使用CentOS6.7单机实验时，对于60G的数据量GeoHash编码完成大概需要18小时，为了提供处理速度，在一个实施例中，采用基于集群的并行化处理方式，基于GeoHash编码的并行化轨迹化简的算法流程参见图7所示，其输入数据是经过预处理之后获得的数据，具体过程是：

第一步、数据分割

将文件系统中(例如HDFS文件系统中的序列化文件)的m个数据文件(每个文件对应一个船舶标识v)划分成m'个数据块Split，每个块交由一个数据节点处理，其中，m'》k，k表示数据节点数。

第二步、Map阶段：并行提取数据

在此步骤中，对数据进行并行GeoHash编码，抽取每条数据的经度x和纬度y并进行GeoHash编码，然后以编码code(编码值)为键，1为键值输出，其形式为<code，1>。

第三步、Reduce阶段：统计不同GeoHash编码code的密度dsy，计算其对应区域C的中心点经纬度C_x和C_y。

经过Reduce阶段的处理之后，以C_x、C_y和dsy为键，null为键值输出，其形式为<(C_x，C_y，dsy)，null>。

第四步、保存数据

所有数据处理完之后，将所有数据按元组<C_x，C_y，dsy>保存。

通过上述处理过程，可以统计具有相同Geohash编码的轨迹点密度(或称轨迹点个数)，例如，以元组<Center.lon，Center.lat，Dsy>保存，其中Center.lon和Center.lat是Geohash编码对应网格的中心点的经度和纬度，Dsy是对应网格内轨迹点密度，图8示出了东海区域栅格化效果示意图。

轨迹化简之后的数据在密度上的分布保持原有的分布特性，为了使相邻栅格密度更加平滑均匀，进一步地，基于图像处理的模糊算法思想，可选择使用修改后的中值滤波算法对结果进行滤波处理，以去除孤立的噪声点并使所有区域密度变化趋势更平滑均匀。

步骤S130，根据栅格内轨迹点的密集程度将较稀疏的栅格进行合并。

在轨迹化简之后，根据轨迹数据分布在一个栅格内的密集程度，对轨迹数据分布比较稀疏的栅格进行合并。

在一个实施例中，采用四叉树进行栅格合并，建立一棵四叉树存储数据，除根节点不存储信息外，每个节点存储一个栅格的GeoHash编码和密度(dsy)两个属性。参见图9示意的四叉树结构，四叉树的建立原理是：从整个世界范围开始，未划分的范围作为根节点(root)，编码为空，然后将整个世界范围分为四块，即经纬度方向各划分一次，得到第一层的四个节点，编码顺序按照00、01、10、11分别对应于西南、东南、西北、东北四个方位，其中，某栅格的编码称为该栅格的GeoHash编码，GeoHash编码的位数简称为“栅格化位数”，一个区域内所划分栅格的栅格化位数的多少反映了该区域的栅格精度；第二层的建立是对第一层的每一块范围进行一分为四，得到第二层的16个节点；类似地，可建立后续层次。通常，所建立的四叉树层次不会超过20层，因为第20层的栅格误差范围是在8米以内，基本满足大多数问题要求的精确度。

在一个优选实施例中，采用精度自适应的合并算法对栅格进行合并，其核心思想是：一定区域内整体栅格密度越高的栅格保留的栅格精度(即栅格化位数)越高，相反，栅格密度越低的栅格保留的栅格精度越低。具体地，首先，设定合并结果的最高栅格精度对应的栅格化位数bitnummax和最低栅格精度对应的栅格化位数bitnummin，以及合并密度阈值dsym；然后，将整个地理范围按照设定的最高栅格精度划分为相等大小的栅格，每个栅格内包含不同个数的轨迹点，对应于栅格的栅格密度值；判断合并过程以四个同属于一个父栅格的子栅格为单位，若这四个子栅格的栅格密度均低于设定的合并密度阈值，那么将该四个子栅格合四为一，即修改它们共同的父栅格的栅格密度为四个子栅格的栅格密度的和，否则，如果不满足合并条件则不进行合并。这种方式采用层次遍历的思想，从最高精度层栅格开始判断，该层栅格合并遍历完成后，再对次高精度层次栅格进行合并操作，对四叉树逐层向上进行合并操作，当合并栅格已经达到设定的最低精度栅格所在层次时，完成合并过程。

例如，首先设定合并结果的最高栅格精度对应的栅格化位数(即栅格化位数)bitnummax为16，最低栅格精度bitnummin为11，合并密度阈值dsym为10，首先，遍历所有相邻4个位数为16的栅格，如果这四个栅格的密度值都小于10，则合并这4个栅格，组成一个15位的栅格，该栅格密度值为4个栅格的密度值之和，否则保留这4个位数为16的栅格；接着，遍历所有相邻4个位数为15的栅格，如果这四个栅格的密度值都小于10，则合并这4个栅格，组成一个14位的栅格，密度值为4个栅格的密度值之和，否则保留这4个位数为15的栅格；以此类推，直到遍历所有相邻4个位数12的栅格，此时得到的栅格最大位数为11；执行完此操作后，所有栅格的位数在16到11之间。如图10所示，在四叉树的实现中，合并密度阈值dsym为10，4个父节点相同的子节点(相邻的四个栅格)的密度分别为dsy0＝3，dsy1＝5，dsy2＝2，dsy3＝4，都小于合并密度阈值，所以合并这四个子节点，合并后父节点的密度为dsy＝dsy0+dsy1+dsy2+dsy3d＝14。

在此实施例中，使用四叉树结构存储数据来进行栅格合并，具有如下优点：该数据结构结合了地理编码划分过程的思想，能够良好的反应数据的层次关系；由于层次一般不超过20层，四叉树结构避免了前缀树的局部遍历查找，所以查找效率高效；通过节点的编码属性得到其上层节点编码，可快速搜索访问到上层节点。

步骤S140，对合并之后的栅格进行过滤，保留属于航道的栅格集。

在进行栅格合并之后，优选地，可对栅格进行过滤，以去除非航道的栅格，保留属于航道的栅格集。

在一个实施例中，采用局部过滤方法。以四叉树窗口局部过滤算法为例，其主要思想是：设定包含固定个数栅格的矩形(称为“窗口”)，从左上方第一个窗口的中心栅格开始，根据栅格的GeoHash编码可得到窗口内相邻的所有栅格的密度值，然后按照一定的固定值或局部过滤公式进行窗口内栅格的过滤计算。窗口的滑动是以窗口中心栅格为滑动对象，计算其相邻的下一个窗口中心栅格，整体按行结构进行滑动，至遍历完成。窗口中心栅格如果和窗口边界的距离大于距离整体范围边界距离，则计算窗口内的相邻点时会出现越界，为了避免计算的邻点越界问题，这里窗口中心点的选取范围则是以一个窗口边沿整个地理范围边对齐后的窗口中心点形成的范围，作为窗口中心点的遍历范围。

例如，采用基于NiBlack二值过滤思想进行局部过滤，过滤阈值用T表示，如果中心点栅格密度大于阈值则保留，否则舍弃。过滤阈值T的计算公式如下：

T＝avg+alpha×var (1)

其中，avg表示窗口内所有元素的平均值，var表示窗口内所有元素的方差，alpha表示方差修正系数。

在统计窗口内所有元素的数值时，某个栅格可能由其更高精度的4的N次方个子栅格构成，这样的栅格则需要统计所有的实际包含的更高精度的栅格的密度值。对于中心栅格同理，如果中心栅格是最低精度栅格，则对它仅包含的一个自身密度值判断是否大于T；否则，如果中心栅格实际保留的是更高精度的子栅格，则对所有子栅格的密度值遍历进行和阈值T的过滤比较。

如图11所示，有一个3*3(用于表示滑动窗口大小)的窗口，窗口中包含不同精度的栅格，其中，中心栅格(包含数字50、10、20和30的区域)又有四个经度更低的子栅格组成，每个栅格中的数字表示该栅格的密度，需要对窗口内所有栅格的密度做计算来决定中心的所有栅格是否保留。

具体地：

第一步、计算窗口中所有栅格密度的均值avg

avg＝(40+10+30+30+50+30+40+50+10+20+30+20+45+35+50)/15＝33

第二步、计算窗口内所有栅格密度的方差var

Var＝[(40-33)²+(10-33)²+(30-33)²+(30-33)²+(50-33)²+(30-33)²+(40-33)²+(50-33)²+(10-33)²+(20-33)²+(30-33)²+(20-33)²+(45-33)²+(35-33)²+(50-33)²]/15＝170

第三步、用Niblack公式计算过滤阈值T。

例如，方差修正系数alpha为0.05，阈值T为：

T＝33+0.05*170＝41

第四步、比较每个中心栅格的密度是否大于过滤阈值T，如果大于则保留该栅格，否则过滤掉该栅格。图11中有四个中心栅格(分别是数字50的区域、数字10的区域、数字20的区域和数字30的区域)，通过与计算得到的T比较，栅格密度大于T的栅格将被保留(即数字50区域)，其他栅格密度都小于T则被过滤。

在进行局部过滤之后，需要进行滑动窗口遍历，如图12所示，左图的黑色网格表示所有需要过滤的栅格，大小为4*4，右图的灰色网格表示滑动窗口，大小为3*3。如果从左上方第一个点开始作为窗口中心栅格，会出现相邻栅格越界的问题，如图13所示，有5个栅格不在整体范围之内，是越界的相邻栅格。因此，优选地，窗口中心点的选取范围是以一个窗口边沿整个地理范围边对齐后的窗口中心点形成的范围，如图14所示，灰色区域的4个栅格是窗口中心点的选择范围。更具体地，如图15所示，窗口中心点遍历的顺序为(1)->(2)->(3)->(4)。

通过上述过程，得到整体过滤后的栅格化结果后，即将不属于航道的栅格过滤掉，得到所有属于航道的栅格集，参见图16示出了东海区域栅格化数据合并和过滤的效果。

栅格化结果在合并过滤之后，可能仍存在原始采集数据的密度分布不均匀特征，所以可以进一步利用图像处理技术中的模糊处理算法，对栅格化结果反复进行滤波，使用周围的栅格密度值的均值代替中心栅格的密度值，达到去除杂乱点、栅格密度分布平滑的目的。

步骤S150，从过滤后获得的栅格集中提取航道区域边界。

在此步骤中，从获得的航道栅格集中正确的提取出相邻栅格所构成的航道区域边界，以及相邻栅格所在区域内部存在的空洞的边界。

在一个实施例中，采用Delaunay三角构网方法来提取航道区域边界。其基本原理是：给定一组平面点集，按照Delaunay方法，可形成Delaunay三角网。

在Delaunay三角网中，三角形的三个点可以构成唯一的一个外接圆，该外接圆半径circum_r的倒数值称为三角形密度指数，表示为：

alpha＝1/circum_r (2)

由于航道内的三角形的密度指数较大，航道外的三角形的密度指数较小，因此，可设定alpha阈值(用alpha_value表示)，如果三角形alpha>＝alpha_value，为航道三角形，保留该三角形的三条边，否则，为空洞三角形，不保留边。

图17示出了利用Delaunay三角网提取航道边界的一个实施例的流程图。输入的是步骤S140中获得的过滤后的数据，首先，使用Delaunay方法对轨迹点集Points进行Delaunay三角化，得到三角面集合Triangles；然后，遍历Triangles，计算每个三角形的密度指数alpha，如果alpha不小于密度指数阈值，则将该三角形的边加入边集Edges；从边集Edges中提取顶点以进行多边形化，构成多边形集Polys；最后，遍历多边形集合Polys，判断每个多边形的面积是否大于面积阈值maxArea，如果大于，则将该多边形的边界坐标polygon.coords添加到Channels集合。最后得出的Channels集合即航道边界多边形的顶点坐标集合。

图18是某海域的一个三角网图，可以看出，航道内轨迹点的密度比较大，这些点构成的三角形的外接圆半径比较小，即三角形密度指数较大，而航道外的轨迹点比较稀疏，这些点构成的三角形的外接圆半径比较大，即三角形密度指数比较小。

图19示出了三角化示意图，图19(a)中有一些轨迹点，对这些轨迹点进行Delaunay构网三角化，即图19(b)所示，可以看到轨迹点比较密集的区域三角形外接圆半径比较小，即密度指数较大，而轨迹点稀疏的地方三角形外接圆半径较大，即密度指数较小，通过设置密度指数阈值可以区分出航道三角形和非航道三角形，图19(c)是提取的航道三角形。

图20是对图19(c)提取的航道三角形进行多边形化的结果。

图21是多边化和过滤的结果，其中左图是多边化后的结果，仍有许多杂乱的面积比较小的多边形，通过设置面积阈值maxArea过滤后，提取的航道边界更加平滑，效果如右图所示。

图22示出了根据本发明的实施例提取的东海区域的航道边界的效果图。

本领域的技术人员应理解的是，上述实施例中涉及的各种阈值，可根据实际进行航道边界提取的区域范围、数据量规模等设置为不同的值，并且某些步骤是优选的，但并不是必须的，例如，去噪过程和通过设置多边形面积阈值过滤掉杂乱的面积比较小的多边形等。

需要说明的是，虽然上文按照特定顺序描述了各个步骤，但是并不意味着必须按照上述特定顺序来执行各个步骤，实际上，这些步骤中的一些可以并发执行，甚至改变顺序，只要能够实现所需要的功能即可。

本发明可以是系统、方法和/或计算机程序产品。计算机程序产品可以包括计算机可读存储介质，其上载有用于使处理器实现本发明的各个方面的计算机可读程序指令。

以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中的技术改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。

Claims (7)

1.一种基于众源轨迹数据的航道信息提取方法，包括以下步骤：

步骤1：将指定的地理范围按照要求的精度划分为多个栅格，众源轨迹数据的每个轨迹点对应到相应的一个栅格，将位于同一个栅格的轨迹点使用该栅格的指定点来表示，获得化简后的数据；

其中，步骤1进一步包括：

对于所述多个栅格进行GeoHash编码，统计不同GeoHash编码的轨迹点密度，计算其对应栅格的中心点的经度和维度，并保存每个栅格的中心点的经度、维度和轨迹点密度作为化简后的数据；

步骤2：根据所述化简后的数据的轨迹点在栅格内的密集程度，按照设定的阈值对栅格进行合并，获得合并后的栅格集；

其中，步骤2进一步包括：

步骤21：将获得的栅格数据利用四叉树结构存储，除根节点之外，每个节点存储一个栅格的GeoHash编码和轨迹点密度；

步骤22：对于四叉树结构中的父节点对应的四个子节点，在该四个子节点的轨迹点密度均小于设定的第一密度阈值时，合并该四个子节点；以及

步骤23：设定包含固定栅格个数的矩形作为窗口，以窗口边沿着整个地理范围边对齐后的窗口中心点形成的范围作为窗口中心点的滑动范围进行遍历，如果窗口中所有栅格的轨迹点密度大于设定的第二密度阈值，则保留该窗口内的栅格；

步骤3：从所述合并后的栅格集提取航道区域边界信息。

2.根据权利要求1所述的方法，其中，所述第二密度阈值设置为：

T＝avg+alpha1×var

其中，avg表示窗口内所有栅格的轨迹点密度的平均值，var表示窗口内所有栅格的轨迹点密度的方差，alpha1表示方差修正系数。

3.根据权利要求1所述的方法，其中，步骤3包括：

步骤31：利用Delaunay对合并后的数据进行三角化，得到三角面集合；

步骤32：计算三角面集合中每个三角形的密度指数，将密度指数小于指数阈值的三角形的边加入到边集；

步骤33：对边集进行多边形化，获得多边形集；

步骤34：将多边形集中面积大于预定的面积阈值的多边形的边界坐标作为航道边界的顶点坐标集合，以形成多边形的航道边界。

4.根据权利要求3所述的方法，其中，所述三角形的密度指数为：

alpha2＝1/circum_r

其中，circum_r表示三角形的外接圆半径。

5.根据权利要求1所述的方法，其中，在步骤1中，所述众源轨迹数据是经过预处理之后的获得的数据，预处理过程包括：

排序步骤：以时间为顺序对众源轨迹数据进行排序；和/或

采样步骤：基于众源轨迹数据的相邻两个轨迹点的时间间隔与第一时间间隔阈值的比较结果确定需要保留的轨迹点；和/或

过滤步骤：根据众源轨迹数据的每个轨迹点的速度与速度阈值的比较结果确定需要过滤掉的轨迹点；和/或

插值步骤：根据过滤后的相邻轨迹点之间的时间间隔与第二时间间隔阈值的比较结果确定在相邻轨迹点之间需要插入的轨迹点。

6.根据权利要求5所述的方法，其中，通过以下步骤执行所述预处理过程：

数据分割步骤：将原始众源轨迹数据划分成m个数据块，每个数据块交由一个数据节点处理，其中m＞＞n，n表示数据节点数；

Map阶段：对于所分配的数据块，多个数据节点并行逐行读取字段无缺失数据，提取每条数据的v、x、y、t四个属性，以字段v为键，元组(x，y，t)为键值输出，其中，v是船舶的唯一标识，x是轨迹点的经度，y是轨迹点的维度，t是时间戳；

Reduce阶段：每个Reduce处理具有相同键v的数据，执行所述排序步骤、采样步骤、过滤步骤和插值步骤。

7.一种计算机可读存储介质，其上存储有计算机程序，其中，该程序被处理器执行时实现根据权利要求1至6中任一项所述方法的步骤。

Images