CN110580388A - 一种基于众源轨迹数据的航道网络提取方法 - Google Patents
一种基于众源轨迹数据的航道网络提取方法 Download PDFInfo
- Publication number
- CN110580388A CN110580388A CN201910801751.6A CN201910801751A CN110580388A CN 110580388 A CN110580388 A CN 110580388A CN 201910801751 A CN201910801751 A CN 201910801751A CN 110580388 A CN110580388 A CN 110580388A
- Authority
- CN
- China
- Prior art keywords
- channel
- triangle
- polygon
- type
- triangles
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000605 extraction Methods 0.000 title claims abstract description 36
- 238000000034 method Methods 0.000 claims abstract description 65
- 238000012545 processing Methods 0.000 claims abstract description 52
- 238000009499 grossing Methods 0.000 claims abstract description 34
- 238000001914 filtration Methods 0.000 claims abstract description 13
- 230000005484 gravity Effects 0.000 claims description 7
- 238000004590 computer program Methods 0.000 claims description 4
- 239000011800 void material Substances 0.000 claims description 4
- 230000008569 process Effects 0.000 claims description 3
- 238000004422 calculation algorithm Methods 0.000 abstract description 14
- 239000000284 extract Substances 0.000 abstract description 8
- 238000010276 construction Methods 0.000 abstract 1
- 230000000694 effects Effects 0.000 description 25
- 238000010586 diagram Methods 0.000 description 16
- 238000012217 deletion Methods 0.000 description 13
- 230000037430 deletion Effects 0.000 description 11
- 239000013598 vector Substances 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 4
- 235000019592 roughness Nutrition 0.000 description 3
- 239000003643 water by type Substances 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 230000033001 locomotion Effects 0.000 description 2
- 238000012544 monitoring process Methods 0.000 description 2
- 230000000717 retained effect Effects 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 230000002159 abnormal effect Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000033228 biological regulation Effects 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 238000013075 data extraction Methods 0.000 description 1
- 238000007418 data mining Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 239000012634 fragment Substances 0.000 description 1
- 239000000446 fuel Substances 0.000 description 1
- 230000006870 function Effects 0.000 description 1
- 238000003780 insertion Methods 0.000 description 1
- 230000037431 insertion Effects 0.000 description 1
- 230000014759 maintenance of location Effects 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 230000001105 regulatory effect Effects 0.000 description 1
- 239000004065 semiconductor Substances 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
- G06Q50/40—Business processes related to the transportation industry
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Business, Economics & Management (AREA)
- Theoretical Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Software Systems (AREA)
- Geometry (AREA)
- Computer Graphics (AREA)
- Economics (AREA)
- General Health & Medical Sciences (AREA)
- Human Resources & Organizations (AREA)
- Marketing (AREA)
- Primary Health Care (AREA)
- Strategic Management (AREA)
- Tourism & Hospitality (AREA)
- General Business, Economics & Management (AREA)
- Processing Or Creating Images (AREA)
Abstract
本发明公开了一种基于众源轨迹数据的航道网络提取方法。本发明通过从AIS中获取大量的船舶轨迹数据中提取航道、对航道边界进行平滑处理、对提取的航道进行三角化及三角形分类处理以过滤掉对航道无用的三角形、对航道进行中心线的提取、航道网络构建的方式实现更加精确的航道网络构建。与现有技术相比，本发明不仅提取海洋航道边界，而且提取中心线和交叉点。本发明提出的方法通过曲率过滤、线性插值、SMA边界平滑算法等各种方式解决了从存在噪声和密度不均匀的轨迹数据中提取完整、连续、平滑的航道网络的问题。
Description
技术领域
本发明涉及数据挖掘技术领域，具体来说，涉及一种基于众源轨迹数据的数据提取方法，更具体地说，涉及一种基于众源轨迹数据的航道网络提取方法。
背景技术
目前，获取地理道路网络信息的传统方式是由训练有素的专职人员在街道上驾驶装有特殊配置的车辆，通过分析高分辨率卫星图像或通过人工众包方式来实现，例如OpenStreetMap就通过众多志愿者的贡献来编辑、完善地图。物联网和云计算技术的进步，使得收集更多的空间轨迹更加方便，例如收集表示汽车、自行车和移动电话等位置信息的数据，这种前所未有的大规模轨迹数据可用于交通监测和预警、轨迹预测、异常检测和旅行路线规划等，在基于位置的服务、城市计算和智能交通系统等领域具有广泛的应用。与传统方法相比，从轨迹数据中提取道路网络更具成本效益、灵活性和更高的时间性能。
在城市计算中，轨迹数据通常从全球定位系统(GPS)终端收集，而在海上交通中，自动识别系统(AIS)终端可以为我们提供可靠的船舶轨迹信息源。虽然海上没有人工建造的实物形式的航道，但船舶出于节省燃油、遵从政府管制等一系列原因，通常在事实上表现为遵循一定的海上航线。因此，如同从GPS数据中提取路网数据一样，人们也有可能从AIS数据中提取出航道。
与基于GPS轨迹数据的航道网络提取相比，从AIS轨迹数据中提取航道网络面临着一系列的挑战，例如：1)AIS数据的数据量大、噪音高、密度和质量不均匀，为航道网络提取提出巨大挑战。长期从大量船只收集的轨迹数据量非常大，并且往往会因为AIS设备采集到的位置数据有时不准确的现象而造成数据质量不好。在我们的全球AIS数据集中，每月有超过100GB的数十亿个轨迹点，超过370000个活跃船舶，每年产生的轨迹点超过1TB。在近海和近岸水域，轨迹点的采样频率在为5秒至100秒和2分钟至10分钟。近海、近岸地区和远海海域AIS数据的密度和质量差异很大。近海和近岸水域的船舶轨迹点比开阔海域更密集。此外，开阔海域的提取精度要求较低，而近海和近岸水域的提取精度要求较高。这些特点对航道网络提取算法提出了巨大挑战。2)面向大范围区域提取连续、平滑的航道网络面临较大挑战。航道网络不仅包括中心线，还需包括边界和交叉点信息等。例如，海洋交通监管系统在监测船只是否航行于合法区域的应用场景中，航道的矢量边界信息具有重要的应用价值。航道的提取方法必须能够保证提取出完整、连续、平滑的航道信息，只能够提取航道片段信息的航道提取方法并不能够在实际中真正应用。传统的聚类算法提取出的航道信息往往只是片段信息，并且应用于大范围海洋轨迹数据时，难以提取出连续、平滑、完整的航道。
现有技术中，从轨迹数据中提取航道网络的相关技术可分为以下三类：(1)基于向量的聚类或几何计算方法：将轨迹点转换为向量，然后提取航道网络的信息。大多数基于向量的方法采用聚类算法进行航道提取。其中，点聚类算法将轨迹点聚集成航道上的关键点，然后将其连接起来形成航道；分段聚类方法通过聚类轨迹段来提取航道。除此之外还有一些基于向量的方法采用几何计算途径，例如，使用移动属性的统计方法来找到航道上的点和线，或通过轨迹的几何关系和/或形状相似性逐步插入以及合并轨迹来提取航道上的线段。上述方法可用于提取航道线或段，但很难用于在大区域内提取航道边界或整个航道边界。此外，该方法适用于密集采样的轨迹数据，对于噪声点多和采样稀疏的轨迹数据效果不好。(2)基于图像的方法：将轨迹点的位置数据转化为图像的横纵坐标，然后使用图像处理技术来识别和/或来提取航道信息。基于图像的方法适用于密集采样轨迹，但对于噪声多的轨迹数据不具有鲁棒性，并且在应用于大量轨迹数据时会有很多问题。(3)基于栅格的方法：将地图转化为若干个同等大小的栅格，计算每个栅格内轨迹点的数量，然后从栅格数据中提取航道信息。有的方法利用Delaunay三角网来提取航道边界和航道中心线，而其他方法使用栅格中的轨迹点统计值来确定道路方向。基于栅格的方法可用于提取边界信息和中心线，并且对于具有噪声和变化密度的数据而言比其他两个方法鲁棒性好，但这些方法未能将其应用于解决基于海量AIS数据、如何在数据质量低且轨迹点密度不均匀的条件下，提取大范围海域完整、连续、平滑的航道。并且，也均没有考虑到航道边界和航道的交叉点的提取。
发明内容
因此，本发明的目的在于克服上述现有技术的缺陷，提供一种新的基于众源轨迹数据的航道网络提取方法。
本发明的一种基于众源轨迹数据的航道网络提取方法，包括如下步骤：
S1、根据船舶的轨迹数据来提取航道；
S2、对步骤S1中提取的航道进行三角化及三角形分类处理；
S3、基于步骤S2处理后的航道数据，进行航道中心线的提取；
S4、根据步骤S3提取的航道中心线，提取航道端结点和航道路段信息，并依据航道端结点和航道路段信息生成航道端结点之间的拓扑关系，再综合航道中心线信息进行航道网络构建。
其中，所述步骤S2包括如下步骤：
S21、对航道进行三角化处理；
S22、对步骤S21中划分的三角形进行过滤，保留三角形重心在航道内部且不在空洞多边形内的三角形；
S23、将经过步骤S22过滤后保留的三角形进行分类，分成4类三角形，分别为0类三角形、1类三角形、2类三角形和3类三角形，其中，0类三角形包含0个公共边、1类三角形包含1个公共边、2类三角形包含2个公共边、3类三角形包含3个公共边；
S24、基于步骤S23中的三角形分类数据进行航道数据处理，将航道中所有包含三角形个数少于10个的以1类三角形开始、以3类三角形结束的航道上的1类三角形和2类三角形删除，并将这些航道上的3类三角形改为2类三角形；
S25、将经过步骤S24处理后的三角形数据进行分类输出，分别输出三角形每条边的信息、三角形每条边的顶点信息、三角形信息及其类型。
经过步骤S25输出的三角形每条边的信息包含该条边两个顶点的标识符以及边的类型、三角形每条边的顶点信息包含三角形每条边对应的两个顶点的经纬度坐标、三角形信息及其类型包含三角形三个顶点的标识符以及三角形的类型。
所述步骤S3包括如下步骤：
S31、基于处理后的航道数据，将航道进行分类，分为第一类航道和第二类航道，第一类航道是以1类三角形开始、以1类或者3类三角形结束的航道，第二类航道是以3类三角形开始、以三类三角形结束的航道；
S32、对分类后的航道进行中心线上的结点提取，得到每条航道中心线上的结点经纬度坐标和组成每条航道的三角形信息，组成航道的三角形具有三角形标识符。
航道中心线上的结点按照以下方式提取：
对于1类三角形，保存其较长的私有边的中点和公共边的中点作为其所在航道中心线的结点；
对于2类三角形，保存其两个公共边的中点作为其所在航道中心线的结点；
对于3类三角形，保存其重心作为其所在航道中心线的结点。
所述步骤S4包括如下步骤：
S41、根据航道中心线的结点经纬度坐标，得到该条航道的航道中心线两个端结点经纬度坐标；
S42、根据组成航道的三角形信息，得到该条航道的中心线两个端结点对应的三角形的类型和三角形标识符；
S43、根据步骤S41获得的航道中心线两个端结点经纬度坐标和步骤S42中获得的航道中心线两个端结点对应的三角形类型和三角形标识符，组合生成航道端结点信息；
S44、根据组成航道的三角形信息，获取航道两端三角形的标识符，由此生成航道路段信息；
S45、根据步骤S43中的航道端结点信息、步骤S44中的航道路段信息，生成航道端结点之间的拓扑关系。
航道端结点之间的拓扑关系是通过如下步骤生成的：
首先，遍历所有航道端结点信息，得到每个端结点对应的三角形标识符信息；
其次，以一个端结点对应的三角形标识符为对象，遍历所有航道路段信息，查看是否有航道路段包含有该三角形标识符，若有，则将该航道路段的另一个端结点对应三角形的标识符加入到邻接表中，以此类推，遍历所有端结点对应的三角形标识符。
提取航道后，对锯齿明显的航道边界进行平滑处理再进行三角化处理。提取的航道是由若干个多边形组成的，包括道路多边形和空洞多边形，航道边界平滑处理包括如下步骤：
H1、提取航道拐点，针对组成航道的每个多边形，通过计算航道多边形的每个点的曲率，将曲率大于第一曲率阈值的点作为拐点；
H2、针对经过步骤H1处理之后的每个航道多边形，计算其短边距离阈值和每个点的曲率阈值，删除与前一个相邻点的距离小于短边距离阈值且曲率大于第二曲率阈值的点；
H3、针对经过步骤H1和H2处理后的航道边界数据进行线性均匀插值处理；
H4、将经过步骤H3插值处理后的航道边界数据采用SMA方式进行平滑处理；
H5、针对航道内的空洞多边形进行凸化处理。
所述步骤H4包括如下步骤：
H41、以组成航道多边形的点的个数为上限，将SMA平滑窗口的大小依次取值为1、3、5、……直到上限，分别计算不同窗口下航道多边形的平滑度和峰度，将航道多边形平滑后峰度大于平滑前峰度且平滑度最小时对应的窗口大小作为航道SMA平滑处理的窗口大小；
H42、采用步骤H41中确定的窗口大小，以组成航道多边形的点作为窗口的中心，计算经纬度的平均值来代替当前点；
H43、针对组成航道多边形的每个点，采用步骤H42进行处理。
所述步骤H5包括如下步骤：
H51、获取航道内的空洞多边形的每个点的经纬度坐标；
H52、针对空洞多边形的每一个点pi，计算该点与其相邻两个点经纬度坐标按照如下方式计算的pi.mul，pi＝(loni，lati)
pi.mul＝(pi-1.loni-1-pi+1.loni+1)(pi+1.lati+1-pi.lati)-
(pi+1.loni+1-pi.loni)(pi-1.lati-1-pi.lati)
若pi.mul大于0则pi-1、pi、pi+1三个点为逆时针排列，将pi点保存在逆时针库中；若pi.mul小于0则相邻三个点为顺时针排列，将pi点保存在顺时针库中；最后对比顺时针库和逆时针库的点的个数，将点的个数多的保留；
H53、针对航道内的每一个空洞多边形，重复步骤H51-H52；
H54、对经过步骤H53处理的空洞多边形进行线性插值处理。
本发明对现有方法进行了改进，提供精确的航道网络提取方法。与现有技术相比，本发明不仅提取海洋航道边界，而且提取中心线和交叉点。本发明提出的方法通过曲率过滤、线性插值、SMA边界平滑算法等各种方式解决了从存在噪声和密度不均匀的轨迹数据中提取完整、连续、平滑的航道网络的问题。
附图说明
以下参照附图对本发明实施例作进一步说明，其中：
图1为根据本发明实施例的一种基于众源轨迹数据的航道网络提取方法流程示意图；
图2为根据本发明实施例采用本发明方法对2016年6月份到12月份的数据提取出的渤海区域的航道示意图；
图3为根据本发明实施例采用本发明方法对图2中的渤海区域的部分航道提取拐点后的效果示意图；
图4为根据本发明实施例采用本发明方法对图2中的渤海区域的部分航道多边形先提取拐点，然后再进行删除曲率较大的短边后的效果对比示意图；
图5为根据本发明实施例采用本发明方法对图2中的渤海区域的部分航道多边形，将提取拐点和短边删除后的数据经过插值处理前后的数据对比示例示意图；
图6为根据本发明实施例采用本发明方法对图2中的渤海区域的部分航道经过提取拐点和删除曲率大的短边后的效果示意图；
图7为根据本发明实施例采用本发明方法对图6中的渤海区域的部分航道进行线性均匀插值的效果示意图；
图8为根据本发明实施例采用本发明方法对图2中的渤海区域的部分航道数据经过平滑后的效果示意图；
图9为根据本发明实施例采用本发明方法对图2中的渤海区域的部分航道的空洞多边形的凸化前后的对比示意图；
图10为根据本发明实施例采用本发明方法对图2中的渤海区域的航道经过航道边界平滑处理后的效果示意图；
图11为根据本发明实施例采用本发明方法对图2中的渤海区域的部分航道经过Delaunay三角化后的效果示意图；
图12为根据本发明实施例采用本发明方法对图11中三角化后的渤海区域的部分航道，将航道外的三角形和空洞内的三角形剔除后的效果示意图；
图13为根据本发明实施例对图2中的渤海区域的部分航道三角化及三角形分类处理后三角形中每一条边的两个顶点的标识符信息；
图14为根据本发明实施例对图2中的渤海区域的部分航道三角化及三角形分类处理后三角形中每一条边的两个顶点的坐标信息；
图15为根据本发明实施例对图2中的渤海区域的部分航道三角化及三角形分类处理后的三角形及其类型的信息；
图16为根据本发明实施例对图2中的渤海区域的部分航道三角化以及三角形分类前的效果示意图；
图17为根据本发明实施例对图2中的渤海区域的部分航道三角化以及三角形分类后的效果示意图；
图18为根据本发明实施例对图2中的渤海区域的航道中心线上结点的经纬度坐标；
图19为根据本发明实施例对图2中的渤海区域的每条航道的三角形的信息；
图20为根据本发明实施例对图2中的渤海区域的航道提取中心线后的效果图；
图21为根据本发明实施例对图2中的渤海区域的航道中心线提取出的航道端结点之间的拓扑关系；
图22为根据本发明实施例对图2中的渤海区域的航道进行航道、航道中心线、航道端结点叠加后的效果示意图。
具体实施方式
为了使本发明的目的，技术方案及优点更加清楚明白，以下结合附图通过具体实施例对本发明进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
为了便于理解，首先介绍本发明涉及的一些概念：航道，是由航道多边形组成的，将航道多边形定义为一个二维多边形形状，代表海上船舶在规范下能够航行的区域范围；道路多边形，即航道的最外侧的多边形；空洞多边形，是由于航道内部存在礁石等障碍物，而出现的不能被船舶行驶的区域范围；非航道，除航道和空洞以外的船舶不允许行驶的区域范围。结点：组成航道中心线上面的点。航道端结点(node)：航道的交叉口或断头路的终点，主要信息包括结点对应三角形的标识符、结点经纬度坐标、结点对应三角形的类型等。航道路段(segment)，是由两个航道端结点组成的边，包括航道路段标识符和起终点端结点对应的三角形标识符。
根据本发明，概括来说，首先从AIS中获取大量的船舶轨迹数据，从轨迹数据中提取航道；然后，对提取的航道进行三角化及三角形分类处理，过滤掉对航道无用的三角形，如果提取的航道锯齿化明显，在三角化之前要进行航道边界平滑处理；第三，对经过三角化和三角形分类处理后的航道进行中心线的提取；最后，根据前面的处理结果，进行航道网络构建。
根据本发明的一个实施例，如图1所示，一种基于众源轨迹数据的航道网络提取方法，包括如下步骤：
根据本发明的一个实施例，本发明的一种基于众源轨迹数据的航道网络提取方法包括如下步骤：
第一步、基于众源轨迹数据的航道提取。根据本发明的一个实施例，航道提取包括以下步骤：对众源轨迹数据的预处理、通过栅格化对预处理后的数据进行化简，获得化简后的轨迹点、根据栅格内轨迹点的密集程度将较稀疏的栅格进行合并、利用Delaunay对合并后的数据进行三角化，得到三角面集合、计算三角面集合中每个三角形的密度指数，将密度指数小于指数阈值的三角形的边加入到边集、对边集进行多边形化，获得多边形集、将多边形集中面积大于预定的面积阈值的多边形的边界坐标作为航道边界的顶点坐标集合，以形成多边形的航道边界；如图2所示，以2016年渤海区域的众源轨迹数据提取的航道结果，提取出的航道是由多个航道多边形组成的，它包含道路多边形和空洞多边形，其中，航道边界有很多锯齿，需要进行平滑处理；
第二步、航道边界平滑处理。由于提取的航道是阶段产生的航道数据，其存在严重的锯齿化问题，所以需要进行平滑处理。航道边界平滑处理分为道路多边形平滑处理和空洞多边形平滑处理，主要是利用SMA算法来达到使航道边界平滑的目的。航道边界平滑处理的过程包括对航道多边形的提取拐点、删除短边、插值和SMA平滑的操作以及对空洞多边形凸化处理的操作，下面根据本发明的一个实施例，详细介绍航道边界平滑处理过程的各个步骤；
1)提取拐点
针对组成航道的航道多边形，计算航道多边形每个点的曲率，对于每个航道多边形，需要把曲率比较小的点进行过滤，设置第一曲率阈值T1，保留大于第一曲率阈值T1的点作为拐点。
我们利用等式1来确定航道多边形当前的点是否是航道边界上的拐点。假设pi-1，pi，pi+1是航道多边形上的三个相邻的点，si表示点pi的曲率。si可以计算为：
在这个等式中，dis(pi,pj)表示点pi和点pj之间的距离。si的值大于等于1；当s的值为1时，此时三点在同一条直线上。si的值越大，pi点的弯曲程度越大；若si>T1，则pi是拐点并将其保存。
其中第一曲率阈值T1的计算方法为T1＝1+beta1，beta1代表着调节因子，取值范围为[0,1]，如果需要保留较多的航道细节，则需将beta值设置的较小。
根据本发明的一个实施例，现有一个航道多边形的一部分点的原始序列，由{p1,p2,p3,...,p15}组成，计算每个点与其前一个相邻点的距离、其后一个相邻点之间距离以及其曲率值，如下表1所示。
在这里，将第一曲率阈值设置为1.031，若点的曲率值小于1.031，则将此点删除，最后得到的拐点序列为{p1,p2,p3,p4,p6,p8,p9,p10,p11,p12}。
表1航道多边形各点间的距离以及各点的曲率值
根据本发明的一个示例，以2016年6月份到12月份的渤海区域的航道为对象，通过提取拐点的方法来过滤掉曲率比第一曲率阈值小的点，在这里第一曲率阈值设置为1.031。渤海区域的部分区域提取拐点前后的对比效果图如图3所示，其中图中折线较多的线代表提取拐点前，直线较多的线代表提取拐点后。
2)删除曲率大的短边
计算每个航道多边形的一个短边距离阈值，并设置第二曲率阈值T2，保留当前点与其相邻的前一个点组成边的长度小于短边距离阈值且其曲率大于第二曲率阈值的点T2，通过删除点来实现短边的删除。其中，第二曲率阈值T2与中第一曲率阈值的计算方法相同,T2＝1+beta2。
经过上述提取拐点处理后，航道多边形现在由拐点组成。由于锯齿状边缘的长度一般比普通边缘短，且当其曲率比较大时，造成航道更加不平滑，我们可以通过删除曲率比较大的短边来使边界更平滑。因此，我们计算航道多边形的每条边的长度，如果一条边的长度小于短边距离阈值并且在此点的曲率大于第二曲率阈值T2，删除这条边。对于每个航道多边形，都需要把曲率大的短边删除。
首先，我们计算每个航道多边形的短边距离阈值ThresDis并设置第二曲率阈值T2。由等式1可确定航道多边形上每个点的曲率值，由等式2可确定每个航道多边形的短边距离阈值ThresDis。其中Mean表示每个航道多边形的边的长度的平均值，stdDev表示边的长度标准差，α表示标准差修正系数并且其取值范围为[0,1]。
ThresDis＝Mean+α*StdDev (2)
假设pi-1，pi，pi+1是航道多边形上的三个相邻的点，计算当前点与其相邻的前一个点之间的距离distancei＝dis(pi,pi-1)和在pi点的曲率si。若distancei小于ThresDis且si>T2，则将点pi删除，以此类推，针对航道多边形上的点按照相同的方式处理。
根据本发明的一个示例，通过对表1进行提取拐点处理，得出的拐点序列{p1,p2,p3,p4,p6,p8,p9,p10,p11,p12}的结果进行短边的删除，计算其每个点与相邻的前一点的距离以及其曲率值，结果如表2所示。
表2：当前点与相邻前一个点的距离以及其曲率值
| 轨迹点编号 | 当前点与前一个相邻点的距离 | 当前点的曲率值 |
| p<sub>1</sub> | -- | -- |
| p<sub>2</sub> | 775.8863 | 1.1839064895162592 |
| p<sub>3</sub> | 1033.9126 | 1.127316596591655 |
| p<sub>4</sub> | 1034.0541 | 1.0170397371196471 |
| p<sub>6</sub> | 1883.2241 | 1.1964810040634752 |
| p<sub>8</sub> | 1601.3755 | 1.3596051749919917 |
| p<sub>9</sub> | 477.3695 | 1.1053335604651593 |
| p<sub>10</sub> | 775.6944 | 1.0460017263849162 |
| p<sub>11</sub> | 305.6746 | 1.0874958148311278 |
| p<sub>12</sub> | 1133.1592 | 1.2153792683561062 |
通过计算得到了短边距离阈值为1500，第二曲率阈值T2设置为1.03。若当前点与相邻前一个点之间的距离小于1500且其曲率值大于1.03，则将其删除，最终得到了以下序列{p1,p4,p6,p8}。
根据本发明的一个实施例，以2016年6月份到12月份的渤海区域的航道为对象，计算出短边距离阈值ThresDis为1500，第二曲率阈值T2为1.03，计算每个点与其相邻的前面一个点的距离以及其曲率，删除距离比1500小并且曲率比1.03大的点，如图4所示，渤海部分区域的删除短边前后的对比效果图如图4所示，其中锯齿较多的线代表删除短边前，平滑弧形的线代表删除短边后。
3)航道多边形线性插值
航道多边形经过提取拐点和删除短边处理以后，数据变得稀疏，影响三角化，本发明在这里利用了线性均匀插值的方法进行插值操作,即计算航道多边形中所有边的平均长度，以平均长度为步长对多边形进行插值；
根据本发明的一个实施例，计算经过提取拐点和删除短边处理之后的多边形中的所有边的总长度dis_sum，计算多边形中所有边的平均长度dis_average＝dis_sum/N-1，其中N代表经过提取拐点和删除短边处理后构成每个航道多边形的点的总个数。计算当前点与其相邻的后面的点之间的距离dis，如果dis大于dis_average，则计算插入点的个数n＝dis/dis_average+1，然后进行均匀插值；否则不处理。
根据本发明的上述示例，对通过删除短边处理之后得到的结果进行插值处理，由于删除短边后多边形的点变为序列{p1,p4,p6,p8}，这只是一个航道多边形的部分数据，我们无法通过这些数据直接得到平均长度dis_average。通过计算这一个航道多边形的平均长度为2278.63773，将序列{p1,p4,p6,p8}通过插值处理得到结果如图5所示。
根据本发明的一个示例，以2016年6月份到12月份的渤海区域的航道为对象，经过提取拐点和删除短边处理过的2016年6月份到12月份的渤海区域的航道数据，在线性均匀插值操作前后的效果如图6和图7所示。
4)SMA平滑
在经过提取拐点、删除短边、线性插值处理以后，航道锯齿边减少了很多，但是还是不够平滑，MA算法或一些基于MA的改进算法可以让航道更加的平滑。在此综合考虑算法的简单、效率及其有效性，本发明的一个实施例采用了SMA算法。SMA算法是通过不断调整窗口的大小，求每个航道多边形的平滑前和平滑后的平滑度(刻画了曲线的平滑程度)和峰度(刻画了曲线对峰值信息的保留程度)，航道边界应该在尽可能保留原有峰值信息的前提下尽可能平滑，因此，在平滑后峰度大于平滑前峰度的前提下，取平滑度最小时的窗口值。以当前点为窗口的中心，计算经纬度的平均值来代替当前点，最后通过窗口的移动来使航道多边形平滑。
经过线性插值处理之后的航道，若要平滑航道多边形，必须要确定需要平滑的航道多边形P以及SMA平滑窗口的大小win_size，其中win_size的值是由等式3和等式4确定的。其中等式3计算了航道多边形的平滑度Roughness，等式4计算了航道多边形的峰度Kurtosis。
Roughness＝σ(ΔS) (3)
在等式3中，ΔS＝(Δs0,Δs1,...)，其中Δsi＝si+1-si，si在等式1中定义，表示的是航道多边形点的曲率，Roughness表示Δs的标准差。
在等式4中，S代表航道多边形P上的点的一系列曲率si值，μ代表S的均值。
经过提取拐点、删除短边、线性插值处理后组成航道多边形的点的个数为channel_size，win_size依次取值为1,3,5,7,...,channel_size。对于win_size取的每个值，将经过提取拐点、删除短边、线性插值处理后得到的航道多边形P平滑后得到P'，分别计算P和P'的Roughness和Kurtosis的值。取Kurtosis[P']大于Kurtosis[P]时且Roughness最小时的win_size值为窗口大小值。
确定了win_size值以后，对于每个航道多边形，设定一个数组result用于记录平滑处理后的点对应的经纬度坐标，计算经过提取拐点、删除短边、线性插值处理后组成每个航道多边形的点个数channel_size，遍历每个航道多边形，其中i代表当前点的下标，取值范围为[0,channel_size-1]，按照以下三个步骤进行处理：
(1)如果channel_size大于等于win_size，则判断i是否小于等于win_size/2，若i小于等于win_size/2，则遍历航道多边形的第0个到第2*i+1个元素，分别计算经度和纬度的平均值并添加到result中；
若i大于win_size/2，则判断i+(win_size/2)是否小于channel_size：若是，则遍历航道多边形的第i-(win_size/2)个到第i+(win_size/2)个元素，分别计算经度和纬度的平均值并进添加到result中；若否，则计算a＝channel_size-i-1，遍历航道多边形的第i-a个到最后一个元素，分别计算经度和纬度的平均值并添加到result中；
(2)如果channel_size小于win_size：则将航道多边形的每个元素添加到result中；
(3)以数组result中记录的经纬度对应的点作为平滑处理后组成航道多边形的点。
根据本发明的一个示例，以2016年6月份到12月份的渤海区域的航道为对象，经过平滑后的部分效果如图8所示，其中锯齿表示平滑前，弧线表示窗口值win_size＝11时SMA平滑后的效果。
5)空洞多边形凸化
由于航道多边形中有些空洞多边形是凹多边形，影响航道中心线的提取，因此需要对空洞多边形进行处理。对空洞多边形的处理是通过删除来使其凸化的，会导致多边形上的点较稀疏，所以最后又进行了线性均匀插值。由于步骤3)已经介绍过线性均匀插值，在此就不多介绍。
对于每个空洞多边形，多边形上的每个点用其经纬度坐标表示为pi＝(loni,lati)，loni表示点pi的经度坐标，lati表示点pi的纬度坐标。首先得到每个点pi的前一个点pi-1以及后一点pi+1的经纬度坐标，通过计算
pi.mul＝(pi-1.loni-1-pi+1.loni+1)(pi+1.lati+1-pi.lati)-(pi+1.loni+1-pi.loni)(pi-1.lati-1-pi.lati)
若pi.mul大于0则pi-1、pi、pi+1三个点为逆时针排列，将pi点保存在逆时针库result1中；若pi.mul小于0则相邻三个点为顺时针排列，将pi点保存在顺时针库result2中。最后分别计算result1和result2中点的个数，将点的个数多的结果进行保存，然后再对保留的点进行线性插值，这样即可将凹多边形变为凸多边形。
根据本发明的一个示例，假设有一个多边形A由以下5个点组成：p1，p2，p3，p4和p5，当前点为p3，将其表示为p3＝(lon3，lat3)，通过计算p3.mul来判断三个点顺序，
p3.mul＝(p2.lon2-p4.lon4)(p4.lat4-p3.lat3)-(p4.lon4-p3.lon3)(p2.lat2-p3.lat3)，
若p3.mul大于0则p2、p3、p4三个点为逆时针排列，将p3点保存在逆时针库result1中；若p3.mul小于0则相邻三个点为顺时针排列，将p3点保存在顺时针库result2中；参照此方法，依次计算多边形A的其他所有点，最后分别计算result1和result2中点的个数，将点的个数多的结果进行保存，针对保留的多边形的点进行线性插值，这样即可将凹多边形变为凸多边形。
根据本发明的一个示例，以2016年6月份到12月份的渤海区域的航道为对象，对2016年6月份到12月份的渤海区域的部分空洞多边形航道数据经过提取拐点、删除短边、线性插值、SMA平滑处理后得到的数据进行空洞多边形凸化处理前后对比的效果图如图9所示，其中内凹的尖角处是凹多边形凸化前，尖角外侧平滑弧线是凹多边形凸化后。2016年6月份到12月份的渤海区域的航道经过航道边界平滑处理后的效果如图10所示，表示了平滑处理之后的结果，其中多折弧线代表航道边界。
第三步，三角化及三角形分类，在进行航道中心线提取之前需要对航道多边形数据进行三角化及三角形的分类，这一步包括以下三角化、三角形过滤、三角形分类和基于三角形分类数据的过滤。进行此步操作的目的是将航道划分为若干个三角形，并将对于航道中心线提取无用的三角形和影响航道中心线提取的较短的航道剔除，以便后面的中心线的提取。其中三角化的目的是为了便于后面将对于提取航道中心线无用的三角形剔除。
在航道中心线提取之前需要对平滑处理后的航道边界进行Delaunay三角剖分和三角形的过滤。通过Delaunay可以将航道划分为若干个三角形，但是会存在一些对于提取航道中心线无用的三角形，如图11中的(2)、(3)、(4)和(5)中的三角形所示。三角形过滤是通过判断三角形的重心是否在航道内部和是否在空洞多边形内来操作的：若三角形重心在航道内部且不在空洞多边形内，则将三角形保留，否则将其过滤。过滤后的效果如图12所示。为了进一步确定中心线上的结点，需要将三角形分类并将某些包含三角形个数少于10个的航道进行删除。其中，三角形分为了以下4类：0类三角形：包含0个公共边；1类三角形：包含1个公共边；2类三角形：包含2个公共边；3类三角形：包含3个公共边。三角形分类之后，再将所有包含三角形个数少与10的以1类三角形开始，以3类三角形结束的航道上的1类三角形和2类三角形删除，并将3类三角形改为2类三角形。
经过处理之后得到了组成每个航道多边形中各个点的标识符id，并输出三种数据，分别命名为side_label、edge_points_filter和triangle_filter_type。在side_label中保存的是所有的三角形中每一条边的两个顶点的标识符id的信息，如图13所示，共由三列组成：前两列分别是这条边的两个顶点的id标识符，第三列是这条边的类型(类型1表示私有边，类型2表示公共边)；在edge_points_filter中保存的是所有三角形中每条边的两个顶点的坐标且其坐标与side_label中id相对应，如图14所示，共由四列组成：前两列表示这条边一个端点的经纬度，后两列表示这条边另一个端点的经纬度；triangle_filter_type保存的是组成航道的三角形及其类型，如图15所示，由四列组成，前三列是三角形的三个顶点的id标识符，第四列表示该三角形的类型。
根据本发明的一个示例，以2016年6月份到12月份的渤海区域的航道为对象，航道多边形经过三角化和三角形分类的前后效果对比图如图16和图17所示，采用2016年渤海区域部分航道为例进行处理，其中图16是航道三角化及三角形分类前的效果，图17是航道三角形化和三角形分类之后的效果。
第四步，对经过平滑处理、三角化及三角形分类处理后的航道进行中心线提取。
在进行三角化和三角形分类处理的过程中，已经将三角形进行了分类。其中，0类三角形是孤立的三角形，不做考虑；1类三角形只会出现在一条航道的两端；3类三角形只存在于多条航道的交叉路口；一条航道除了两个端点其他都是2类三角形。根据三角形分类，找出第一类航道和第二类航道，然后将组成每条航道的三角形的信息保存到每条航道的三角形列表中trianglelist中，并针对不同类型的三角形做不同的处理，分别提取每条航道中心线上的结点到中心线列表centerlinelist中。其中，第一类航道是指以1类三角形开始，1类或3类三角形结束的航道；第二类航道是以3类三角形开始，3类三角形结束的航道。
对于每条航道中心线上的结点是由以下方式提取的：对于1类三角形，保存其较长私有边的中点和公共边的中点作为其对应航道中心线的一个结点；对于2类三角形，保存其两个公共边的中点作为其对应航道中心线的两个结点；对于3类三角形，保存其重心作为其对应航道中心线的一个结点。
在提取中心线后，得到的数据为centerlinelist和trianglelist，其中，centerlinelist表示的是航道中心线上结点的经纬度坐标，如图18所示；trianglelist表示的是组成每条航道的三角形的信息，如图19所示，共由8列组成，前6列表示组成每条航道的三角形三个顶点的经纬度，第7列是三角形类型，第8列是三角形trid标识符，组成航道的每一个三角形有唯一的标识符trid。
经过本步骤的处理之后，输出每条航道中心线数据和组成每条航道的三角形信息。
其中，每条航道的中心线由组成该中心线的所有结点的来表示，表示为：
centerlinelistx＝{Pointx1,Pointx2,...,Pointxn}
其中centerlinelistx代表标识符为x的航道中心线，Pointxj＝(lonxj,latxj)代表组成标识符为x的航道中心线上的第j个结点，结点由其经纬度坐标表示；
组成每条航道的三角形的信息，表示为：
trianglelistx＝{trianglex1,trianglex2,...,trianglexm}
其中trianglelistx代表组成标识符为x的航道的所有三角形，trianglexk＝(lonxk1,latxk1,lonxk2,latxk2,lonxk3,latxk3,typexk,tridxk)代表组成标识符为x的航道的第k个三角形的三个顶点的经纬度坐标、三角形类型typexk和三角形标识符tridxk。
根据本发明的一个示例，以2016年6月份到12月份的渤海区域的航道为对象，提取渤海区域航道中心线如图20所示。
第五步，航道网络构建。经过第四步的处理得到了航道中心线，但只有中心线不能清晰地表达海洋路网各要素的特征和相互关系，也不能清楚地表达海洋路网各要素的特征和相互关系，并且不易于用于各种应用，如轨迹预测和路径规划。因此，在本步骤中，提取了航道中心线上的航道端结点node、航道路段segment和航道端结点之间的拓扑关系。根据所有航道端结点之间形成的拓扑关系，可以运用于之后的轨迹预测中，当给出船舶当前所在的航道的交叉路口时，可以很方便地查到可选择航行的路段并且进行下一步的航行方向预测与规划。
根据第四步处理得到的航道中心线数据和组成航道的三角形信息，可以得到每个端结点、路段和所有端结点之间的拓扑关系，根据本发明的一个实施例，其包括如下步骤：
1)生成端结点
根据centerlinelistx的Pointx1和Pointxn可以得到标识符为x的航道的中心线的两端的坐标，根据trianglelistx的第一个元素trianglex1的typex1和tridx1和最后一个元素trianglexm的typexm和tridxm，可以得到标识符为x的航道中心线的两端的结点所在三角形的类型和此结点所在三角形的标识符。其中typexk取值只能为1或3，标识符为x的航道端结点表示如下：
nodeMapx1＝(tridx1,lonx1,latx1,typex1)
nodeMapx2＝(tridxm,lonxn,latxn,typexm)
2)生成航道路段
根据trianglelistx的tridx1和tridxm可以得到组成标识符为x的航道的两端的三角形的标识符。标识符为x的航道的航道路段用edgelistx＝(x,tridx1,tridxm)表示；
3)生成端结点之间的拓扑关系
遍历航道中所有航道端结点nodeMap中的所有元素，得到每个航道端结点对应的三角形标识符trid，通过遍历航道路段edgelist来查看是否存在一条航道路段包含此端结点对应三角形的trid；若包含则将此条航道路段的另一个端结点对应的标识符trid加入到邻接表中。其中邻接表的结构如图21所示，第一列表示一个端结点，其他列表示与第一列的端结点连接的端结点。
根据本发明的一个示例，现有5条航道，其标识符序列为{0,1,2,3,4}，从这几条航道中得到了各条航道的三角形信息和其中心线，分别如表3和表4所示：
表3组成每条航道的三角形信息
表4每条航道的中心线信息
根据表3和表4可得到各个航道的端结点和路段的信息，如表5所示。
表5各条航道的端结点和路段信息
| 航道标识符 | 端结点 | 路段 |
| 0 | {(10,lon<sub>01</sub>,lat<sub>01</sub>,1),(12,lon<sub>04</sub>,lat<sub>04</sub>,3)} | (0,10,12) |
| 1 | {(13,lon<sub>11</sub>,lat<sub>11</sub>,1),(12,lon<sub>04</sub>,lat<sub>04</sub>,3)} | (1,13,12) |
| 2 | {(12,lon<sub>04</sub>,lat<sub>04</sub>,3),(18,lon<sub>25</sub>,lat<sub>25</sub>,3)} | (2,12,18) |
| 3 | {(18,lon<sub>25</sub>,lat<sub>25</sub>,3),(20,lon<sub>34</sub>,lat<sub>34</sub>,1)} | (3,18,20) |
| 4 | {(18,lon<sub>25</sub>,lat<sub>25</sub>,3),(22,lon<sub>44</sub>,lat<sub>44</sub>,1)} | (4,18,22) |
根据表5可生成邻接表如表6所示。
表6：由各个端结点生成的邻接表信息
根据本发明的一个示例，以2016年6月份到12月份的渤海区域的航道为对象，根据2016年的6月份到12月份数据，得到了渤海区域的航道的航道端结点、航道路段和航道端结点之间的拓扑关系。端结点之间的拓扑关系如图21所示。将渤海区域的平滑之后的航道、航道中心线和每条航道的端结点三者叠加显示于图22中。
本发明对现有方法进行了改进，提供精确的航道网络提取方法。与现有技术相比，本发明不仅提取海洋航道边界，而且提取中心线和交叉点。本发明提出的方法是使用基于并行计算的算法有效地处理大量轨迹数据，并且通过栅格合并、滑动窗口过滤、SMA边界平滑算法等各种方式解决了从存在噪声和密度不均匀的轨迹数据中提取完整、连续、平滑的航道网络的问题。
需要说明的是，虽然上文按照特定顺序描述了各个步骤，但是并不意味着必须按照上述特定顺序来执行各个步骤，实际上，这些步骤中的一些可以并发执行，甚至改变顺序，只要能够实现所需要的功能即可。
本发明可以是系统、方法和/或计算机程序产品。计算机程序产品可以包括计算机可读存储介质，其上载有用于使处理器实现本发明的各个方面的计算机可读程序指令。
计算机可读存储介质可以是保持和存储由指令执行设备使用的指令的有形设备。计算机可读存储介质例如可以包括但不限于电存储设备、磁存储设备、光存储设备、电磁存储设备、半导体存储设备或者上述的任意合适的组合。计算机可读存储介质的更具体的例子(非穷举的列表)包括：便携式计算机盘、硬盘、随机存取存储器(RAM)、只读存储器(ROM)、可擦式可编程只读存储器(EPROM或闪存)、静态随机存取存储器(SRAM)、便携式压缩盘只读存储器(CD-ROM)、数字多功能盘(DVD)、记忆棒、软盘、机械编码设备、例如其上存储有指令的打孔卡或凹槽内凸起结构、以及上述的任意合适的组合。
以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中的技术改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。
Claims (14)
1.一种基于众源轨迹数据的航道网络提取方法，其特征在于，包括如下步骤：
S1、根据船舶的轨迹数据来提取航道；
S2、对步骤S1中提取的航道进行三角化及三角形分类处理；
S3、基于步骤S2处理后的航道数据，进行航道中心线的提取；
S4、根据步骤S3提取的航道中心线，提取航道端结点和航道路段信息，并依据航道中心线、航道端结点和航道路段信息生成航道端结点之间的拓扑关系，进行航道网络构建。
2.根据权利要求1所述的一种基于众源轨迹数据的航道网络提取方法，其特征在于，
所述步骤S2包括如下步骤：
S21、对航道进行三角化处理；
S22、对步骤S21中划分的三角形进行过滤，保留三角形重心在航道内部且不在空洞多边形内的三角形；
S23、将经过步骤S22过滤后保留的三角形进行分类，分成4类三角形，分别为0类三角形、1类三角形、2类三角形和3类三角形，其中，0类三角形包含0个公共边、1类三角形包含1个公共边、2类三角形包含2个公共边、3类三角形包含3个公共边；
S24、基于步骤S23中的三角形分类数据进行航道数据处理，将航道中所有包含三角形个数少于10个的以1类三角形开始、以3类三角形结束的航道上的1类三角形和2类三角形删除，并将这些航道上的3类三角形改为2类三角形；
S25、将经过步骤S24处理后的三角形数据进行分类输出，分别输出三角形每条边的信息、三角形每条边的顶点信息、三角形信息及其类型。
3.根据权利要求2所述的一种基于众源轨迹数据的航道网络提取方法，其特征在于，经过步骤S25输出的三角形每条边的信息包含该条边两个顶点的标识符以及边的类型、三角形每条边的顶点信息包含三角形每条边对应的两个顶点的经纬度坐标、三角形信息及其类型包含三角形三个顶点的标识符以及三角形的类型。
4.根据权利要求3所述的一种基于众源轨迹数据的航道网络提取方法，其特征在于，所述步骤S3包括如下步骤：
S31、基于处理后的航道数据，将航道进行分类，分为第一类航道和第二类航道，第一类航道是以1类三角形开始、以1类或者3类三角形结束的航道，第二类航道是以3类三角形开始、以三类三角形结束的航道；
S32、对分类后的航道进行中心线上的结点提取，得到每条航道中心线上的结点经纬度坐标和组成每条航道的三角形信息，组成航道的三角形具有三角形标识符。
5.根据权利要求4所述的一种基于众源轨迹数据的航道网络提取方法，其特征在于，航道中心线上的结点按照以下方式提取：
对于1类三角形，保存其较长的私有边的中点和公共边的中点作为其所在航道中心线的结点；
对于2类三角形，保存其两个公共边的中点作为其所在航道中心线的结点；
对于3类三角形，保存其重心作为其所在航道中心线的结点。
6.根据权利要求5所述的一种基于众源轨迹数据的航道网络提取方法，其特征在于，所述步骤S4包括如下步骤：
S41、根据航道中心线的结点经纬度坐标，得到该条航道的航道中心线两个端结点经纬度坐标；
S42、根据组成航道的三角形信息，得到该条航道的中心线两个端结点对应的三角形的类型和三角形标识符；
S43、根据步骤S41获得的航道中心线两个端结点经纬度坐标和步骤S42中获得的航道中心线两个端结点对应的三角形类型和三角形标识符，组合生成航道端结点信息；
S44、根据组成航道的三角形信息，获取航道两端三角形的标识符，由此生成航道路段信息；
S45、根据步骤S43中的航道端结点信息、步骤S44中的航道路段信息，生成航道端结点之间的拓扑关系。
7.根据权利要求6所述的一种基于众源轨迹数据的航道网络提取方法，其特征在于，航道端结点之间的拓扑关系是通过如下步骤生成的：
首先，遍历所有航道端结点信息，得到每个端结点对应的三角形标识符信息；
其次，以一个端结点对应的三角形标识符为对象，遍历所有航道路段信息，查看是否有航道路段包含有该三角形标识符，若有，则将该航道路段的另一个端结点对应三角形的标识符加入到邻接表中，以此类推，遍历所有端结点对应的三角形标识符。
8.根据权利要求1所述的一种基于众源轨迹数据的航道网络提取方法，其特征在于，提取航道后，对航道边界进行平滑处理再进行三角化处理。
9.根据权利要求8所述的一种基于众源轨迹数据的航道网络提取方法，其特征在于，提取的航道是由若干个多边形组成的，包括道路多边形和空洞多边形，航道边界平滑处理包括如下步骤：
H1、提取航道拐点，针对组成航道的每个多边形，通过计算航道多边形的每个点的曲率，将曲率大于第一曲率阈值的点作为拐点；
H2、针对经过步骤H1处理之后的每个航道多边形，计算其短边距离阈值和每个点的曲率阈值，删除与前一个相邻点的距离小于短边距离阈值且曲率大于第二曲率阈值的点；
H3、针对经过步骤H1和H2处理后的航道边界数据进行线性均匀插值处理；
H4、将经过步骤H3插值处理后的航道边界数据采用SMA方式进行平滑处理；
H5、针对航道内的空洞多边形进行凸化处理。
10.根据权利要求9所述的一种基于众源轨迹数据的航道网络提取方法，其特征在于，所述步骤H4包括如下步骤：
H41、以组成航道多边形的点的个数为上限，将SMA平滑窗口的大小依次取值为1、3、5、……直到上限，分别计算不同窗口下航道多边形的平滑度和峰度，将航道多边形平滑后峰度大于平滑前峰度且平滑度最小时对应的窗口大小作为航道SMA平滑处理的窗口大小；
H42、采用步骤H41中确定的窗口大小，以组成航道多边形的点作为窗口的中心，计算经纬度的平均值来代替当前点；
H43、针对组成航道多边形的每个点，采用步骤H42进行处理。
11.根据权利要求10所述的一种基于众源轨迹数据的航道网络提取方法，其特征在于，所述步骤H5包括如下步骤：
H51、获取航道内的空洞多边形的每个点的经纬度坐标；
H52、针对空洞多边形的每一个点pi，计算该点与其相邻两个点经纬度坐标按照如下方式计算的pi.mul，pi＝(loni，lati)
pi.mul＝(pi-1.loni-1-pi+1.loni+1)(pi+1.lati+1-pi.lati)-
(pi+1.loni+1-pi.loni)(pi-1.lati-1-pi.lati)
若pi.mul大于0则pi-1、pi、pi+1三个点为逆时针排列，将pi点保存在逆时针库中；若pi.mul小于0则相邻三个点为顺时针排列，将pi点保存在顺时针库中；最后对比顺时针库和逆时针库的点的个数，将点的个数多的保留；
H53、针对航道内的每一个空洞多边形，重复步骤H51-H52。
12.根据权利要求9所述的一种基于众源轨迹数据的航道网络提取方法，其特征在于，还包括如下步骤：
H54、对经过步骤H53处理的空洞多边形进行线性插值处理。
13.一种计算机设备，包括存储器和处理器，在所述存储器上存储有能够在处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序时实现权利要求1至12中任一项所述的方法的步骤。
14.一种存储介质，存储有能够在处理器上运行的计算机程序，所述程序运行时执行权利要求1至12中任一项所述的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910801751.6A CN110580388B (zh) | 2019-08-28 | 2019-08-28 | 一种基于众源轨迹数据的航道网络提取方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910801751.6A CN110580388B (zh) | 2019-08-28 | 2019-08-28 | 一种基于众源轨迹数据的航道网络提取方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110580388A true CN110580388A (zh) | 2019-12-17 |
| CN110580388B CN110580388B (zh) | 2023-04-07 |
Family
ID=68812053
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910801751.6A Active CN110580388B (zh) | 2019-08-28 | 2019-08-28 | 一种基于众源轨迹数据的航道网络提取方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110580388B (zh) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112699196A (zh) * | 2020-12-30 | 2021-04-23 | 深圳赛安特技术服务有限公司 | 轨迹生成方法、装置、终端设备及存储介质 |
| CN113312438A (zh) * | 2021-03-09 | 2021-08-27 | 中南大学 | 融合航线提取与趋势判断的海上目标位置预测方法 |
| CN113538482A (zh) * | 2021-06-24 | 2021-10-22 | 电子科技大学长三角研究院(湖州) | 一种基于浮动车轨迹位置和方向信息的车道级城市道路提取方法 |
| CN114019967A (zh) * | 2021-10-29 | 2022-02-08 | 中国船舶重工集团公司第七0七研究所 | 一种适用于狭长航道的无人艇航线规划方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7002578B1 (en) * | 1999-05-05 | 2006-02-21 | Siemens Aktiengellschaft | Method of obtaining a three-dimensional map representation and navigation system |
| US7653218B1 (en) * | 2006-05-02 | 2010-01-26 | Orbimage Si Opco, Inc. | Semi-automatic extraction of linear features from image data |
| US20130328863A1 (en) * | 2012-06-10 | 2013-12-12 | Apple Inc. | Computing plausible road surfaces in 3d from 2d geometry |
| CN106528740A (zh) * | 2016-11-04 | 2017-03-22 | 中科宇图科技股份有限公司 | 基于Delaunay三角网的道路中心线提取方法 |
| CN108920481A (zh) * | 2018-04-20 | 2018-11-30 | 中国地质大学（武汉） | 一种基于手机定位数据的道路网重建方法及系统 |
-
2019
- 2019-08-28 CN CN201910801751.6A patent/CN110580388B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7002578B1 (en) * | 1999-05-05 | 2006-02-21 | Siemens Aktiengellschaft | Method of obtaining a three-dimensional map representation and navigation system |
| US7653218B1 (en) * | 2006-05-02 | 2010-01-26 | Orbimage Si Opco, Inc. | Semi-automatic extraction of linear features from image data |
| US20130328863A1 (en) * | 2012-06-10 | 2013-12-12 | Apple Inc. | Computing plausible road surfaces in 3d from 2d geometry |
| CN106528740A (zh) * | 2016-11-04 | 2017-03-22 | 中科宇图科技股份有限公司 | 基于Delaunay三角网的道路中心线提取方法 |
| CN108920481A (zh) * | 2018-04-20 | 2018-11-30 | 中国地质大学（武汉） | 一种基于手机定位数据的道路网重建方法及系统 |
Non-Patent Citations (3)
| Title |
|---|
| WEI YANG 等: ""A Method for Extracting Road Boundary Information from Crowsourcing Vehicle GPS Trajrctories"", 《SENSORS》 * |
| 李功权 等: ""一种基于约束三角网的道路中心线的提取方法"", 《长江大学学报（自科板）》 * |
| 王振华 等: ""众源GPS浮动车数据中城市道路中心线分级提取的栅格化方法"", 《测绘通报》 * |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112699196A (zh) * | 2020-12-30 | 2021-04-23 | 深圳赛安特技术服务有限公司 | 轨迹生成方法、装置、终端设备及存储介质 |
| CN112699196B (zh) * | 2020-12-30 | 2023-10-27 | 深圳赛安特技术服务有限公司 | 轨迹生成方法、装置、终端设备及存储介质 |
| CN113312438A (zh) * | 2021-03-09 | 2021-08-27 | 中南大学 | 融合航线提取与趋势判断的海上目标位置预测方法 |
| CN113312438B (zh) * | 2021-03-09 | 2023-09-15 | 中南大学 | 融合航线提取与趋势判断的海上目标位置预测方法 |
| CN113538482A (zh) * | 2021-06-24 | 2021-10-22 | 电子科技大学长三角研究院(湖州) | 一种基于浮动车轨迹位置和方向信息的车道级城市道路提取方法 |
| CN113538482B (zh) * | 2021-06-24 | 2024-02-02 | 电子科技大学长三角研究院(湖州) | 一种基于浮动车轨迹位置和方向信息的车道级城市道路提取方法 |
| CN114019967A (zh) * | 2021-10-29 | 2022-02-08 | 中国船舶重工集团公司第七0七研究所 | 一种适用于狭长航道的无人艇航线规划方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110580388B (zh) | 2023-04-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110580388B (zh) | 一种基于众源轨迹数据的航道网络提取方法 | |
| Liu et al. | Adaptive douglas-peucker algorithm with automatic thresholding for AIS-based vessel trajectory compression | |
| CN109241069B (zh) | 一种基于轨迹自适应聚类的路网快速更新的方法及系统 | |
| CN110095127B (zh) | 一种基于分段的隐马尔可夫模型地图匹配方法 | |
| Ai et al. | A DEM generalization by minor valley branch detection and grid filling | |
| CN111291776B (zh) | 基于众源轨迹数据的航道信息提取方法 | |
| Uduwaragoda et al. | Generating lane level road data from vehicle trajectories using kernel density estimation | |
| CN107610223A (zh) | 基于LiDAR点云的电力塔三维重建方法 | |
| CN113450596A (zh) | 基于船舶轨迹特征点提取的时空dp方法 | |
| CN109544691B (zh) | 自动融合多源异构水深数据构建高分辨率dbm的mf方法 | |
| CN111858810A (zh) | 一种面向道路dem构建的建模高程点筛选方法 | |
| CN110929415A (zh) | 一种针对高原地区u型河谷水下未知地形精确模拟方法 | |
| CN111062958A (zh) | 一种城市道路要素提取方法 | |
| CN116543310A (zh) | 一种基于Voronoi图和核密度的道路线提取方法 | |
| CN112116709A (zh) | 一种提高地形表达精度的地形特征线处理方法 | |
| Mangala et al. | A new automatic road extraction technique using gradient operation and skeletal ray formation | |
| CN111080080A (zh) | 一种村镇地质灾害风险预估方法及系统 | |
| CN110411450A (zh) | 一种用于压缩轨迹的地图匹配方法 | |
| Qiu et al. | Inferring road maps from sparsely-sampled GPS traces | |
| Meratnia et al. | A new perspective on trajectory compression techniques | |
| CN110737931A (zh) | 一种基于ArcGIS的铁路桥渡水文关键参数提取方法 | |
| Li et al. | The parallel and precision adaptive method of marine lane extraction based on QuadTree | |
| Sandro et al. | Model generalization of the hydrography network in the CARGEN project | |
| CN114792397A (zh) | 一种sar影像城市道路提取方法、系统以及存储介质 | |
| Liu et al. | Extraction of hottest shipping routes: From positioning data to intelligent surveillance |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |