海运是国际贸易中最主要的运输方式[1]．近年来，随着世界经济的快速发展，各国港口有数量空前的锚泊及排队等待卸货的船舶．船舶锚地和锚位需求的日益增长与锚地资源相对紧缺的矛盾也愈发突出[2]．为了提高船舶锚泊的效率与安全，部分学者从锚位面积、锚地容量、锚位规划等方面出发开展了相关研究．周春辉等[3]提出一种基于Voronoi图的内河船舶锚位面积计算方法，计算和分析锚泊船舶占用水域面积情况．唐皇等[4]采用聚类方法从船舶轨迹数据中提取船舶锚泊数据，并使用关联规则算法提取船舶锚泊半径．郭子坚等[5-6]对锚地实用容量进行概率评价，并对锚地冗余面积进行优化．刘明俊等[7]选取多种代表船型并对其在风、流影响下所需的锚位面积进行计算，最终运用排队论对锚地容量的大小进行动态分析．Oz等[8]提出了一种基于多目标优化策略的船舶锚泊位置选择方法，旨在提高锚地利用率和减少锚泊事故的风险．Madadi等[9]提出了一种考虑锚地时空利用率、船舶碰撞风险的多目标锚地规划算法，仿真实验表明该方法可以优化锚地利用率和降低事故风险．综上，无论是锚地容量计算还是锚位规划等都离不开船舶锚位面积的计算，因此如何合理地计算船舶锚位面积，对于优化锚地资源、提高锚泊的效率和安全至关重要．同时，船舶锚位面积的正确计算对于加强锚泊船舶的安全监管，降低船舶走锚和碰撞的风险有着积极意义．船舶自动识别系统(automatic identification system，AIS)，是指一种应用于船和岸、船和船之间的海事安全与通信的助航系统．通过AIS可以轻松获取船舶的船名、载货、吃水等重要信息[10]．因此，由AIS数据形成的船舶轨迹可以高精度地还原船舶的锚泊情况．本研究首先对AIS轨迹数据进行集成、清洗、分段等预处理，其次使用一种基于密度的船舶停留点识别方法提取船舶锚位信息，再次使用单调链凸包算法求取锚位边界，最后分别使用鞋带定理与多边形形心公式完成船舶锚位面积与落锚点概位的求取．1 锚位面积计算方法1.1　数据预处理低质量的AIS数据将导致低质量的提取结果．为提高数据质量，须要对数据进行预处理．本研究数据预处理主要包括数据集成、数据清洗、轨迹分段．a．数据集成船舶锚泊分析中只有位置报告是不够的，船舶静态和航行相关数据必不可少，因此须要将不同报文数据整合在一起，即数据集成．本研究依据AIS报文共有属性水上移动通信业务标识码(maritime mobile service identify，MMSI)实现不同报文的匹配，集成数据．b．数据清洗AIS数据具有高采样率和异源性，数据中不可避免地会出现冗余和噪声，因此须要对数据进行清洗[11]．AIS数据清洗主要包括去重、光滑噪声等[12]．对于噪声点的发现与处理，本研究选用了一种较为简单有效的处理方式．图1为AIS轨迹去噪示意图，图中：p1~p9为船舶轨迹点，其连线为原始轨迹．为了找出p5这个噪声轨迹点，首先计算点p4到点p5之间的距离Δd，再除以p4和p5两点的时间间隔Δt，得到p4到p5之间的平均速度va．当va大于设定阈值vth时，例如40 kn (全球商船的平均速度不足12 kn，最大速度一般不超过30 kn，当计算速度值超过40 kn时，基本可以认定为异常)，判定p5为异常轨迹点．找到异常轨迹点后，直接删除该噪声点，从而完成轨迹去噪处理．10.13245/j.hust.239467.F001图1AIS轨迹去噪示意图c．轨迹分段图2为AIS轨迹分段流程．当获取一条特定区域、特定时间段的船舶轨迹时，轨迹大多时候是比较混乱的，包含自相交、重复行驶等情况，并不利于做数据分析．因此，须要将一个完整轨迹按照规律切分为多个不同的轨迹段，以便于进一步分析．选取时间间隔作为分段条件对轨迹进行分段．考虑AIS的航行状态报告的报告间隔为2~180 s，一般而言当连续轨迹点间的时间差大于600 s时，即可认定该轨迹点为分段点．10.13245/j.hust.239467.F002图2AIS轨迹分段流程1.2　锚位识别船舶锚位是指一个可供单一船只、海上飞机等停泊的指定水域，一般来说，是由中心点(落锚点)和一个旋回圈来定义．本研究通过AIS轨迹数据完成船舶锚位识别．虽然船舶AIS位置报告具有导航状态信息，但其设定须人工完成．在实际操作过程中，存在参数修改不及时或设定错误等问题[11]．通过导航状态参数获取船舶锚泊信息是不准确的，因此须要对船舶轨迹进行锚位识别，提取船舶锚泊水域范围．本研究以文献[13]为基础，考虑船舶轨迹点时空聚集的同时，兼顾船舶轨迹点的时间连续性和方向性，提出一种基于密度的船舶锚位提取方法．该方法的基本流程如图3所示，主要包括三个步骤：密度计算与候选停留点列表生成，停留点识别，锚位识别．10.13245/j.hust.239467.F003图3船舶锚位提取流程a．密度计算与候选停留点列表生成首先，依次以轨迹点序列中的各轨迹点为锚点，沿时间轴后向搜索和前向搜索，得到各锚点的密度区间(如给定船舶轨迹点序列中存在一最长连续子序列，其中任意一点与该锚点间的球面距离小于设定距离阈值，则该子序列为该锚点的密度区间)．其次，计算各锚点的密度(其密度区间中轨迹点数目)和时间跨度(其密度区间起点与终点的时间差)．再次，判断各锚点的时间跨度是否大于设定时间阈值，是则该锚点为候选停留点．最后，生成候选停留点列表．其中锚点的密度计算原理如图4所示，图中：p1~p15为船舶轨迹点，p8为当前锚点，p8的密度区间M(p8)为p5到p11的连续轨迹子序列．10.13245/j.hust.239467.F004图4密度计算示意图b．停留点识别首先，依据密度对候选停留点列表进行降序排列，选取第一个候选停留点作为当前停留点pstay；其次，遍历候选停留点列表，去除各候选停留点的密度区间与pstay的密度区间的相交部分；再次，将pstay加入到停留点集合，并且重新对列表中各候选停留点进行密度和时间跨度计算；然后，将列表中时间跨度小于设定时间阈值的点进行剪枝，更新列表；最后，重复以上各步，直至候选停留点列表为空．c．锚位识别船舶停留有靠泊和锚泊两种状态．在锚泊过程中，由于水动力、风动力和锚链张力的影响，因此船舶会出现周期性偏荡．特别是当船舶锚泊在潮汐水域时，船舶的位置会随着潮流的变化而变化．因此，停留点位置与速度的标准差在锚泊期间较大，在靠泊期间较小．计算每个停留点位置和速度的标准差，若标准差大于设定阈值，则为锚位点；否则，它是一个靠泊点．1.3　锚位边界提取船舶锚泊时段内的轨迹点是离散的、无明显规则的．通过1.2节计算得来的锚泊点须进行边界提取，以获取船舶实际的锚泊水域范围．对于散点轮廓的提取，目前已有多种方法被提出，比较典型的有凸包(最外层点连接起来构成的一个包含所有给定点的最小凸多边形)与凹包(最外层点连接起来构成的一个包含所有给定点的最小非凸多边形)算法[14]．然而对于同一组散点数据，两种算法可能会得到不同的结果．如图5所示，相较于凹包算法，凸包算法所求得轮廓范围较大．考虑到锚泊船舶运动的随机性与AIS数据的精度，更大的锚泊水域范围，意味着更多的旋回水域，更大的安全余量．因此本研究选取凸包算法求取锚泊点边界，力求通过AIS数据将船舶锚泊水域范围限定到最小，但又最大限度地保留安全余量．10.13245/j.hust.239467.F005图5散点轮廓提取本研究使用时间复杂度为O的单调链凸包(monotone chain convex hull)算法[15]求解锚泊点边界．算法大致流程如下．首先，把所有锚泊点以经度坐标为第一关键字，纬度坐标为第二关键字排序，排序后最小的元素和最大的元素一定在凸包上．然后，在寻找凸多边形过程中，如果从最左边的轨迹点出发顺时针走，轨迹总是右拐的(一旦出现左拐，就说明这一段不在凸包上)，最终到达最右边的轨迹点．这些轨迹点将形成上凸壳．同理从最右边点开始，尝试顺时针方向到达最左边，形成下凸壳．最后，将上下凸壳组合并返回凸包，即锚位的边界．1.4　锚位面积与落锚点计算受风流影响，船舶锚泊期间的活动范围是不规则的．通过锚位边界提取算法求得的锚位边界轮廓是个不规则凸多边形．因此，锚位面积的求取变成为对任意凸多边形面积的求取．本研究使用鞋带公式[16]完成对锚位面积的求取．鞋带公式是在给定顶点坐标的情况下计算平面中非自相交多边形面积的方法．该公式计算过程如图6所示．10.13245/j.hust.239467.F006图6鞋带公式计算过程示意图鞋带公式的表达式为S=12x1y2+x2y3+⋯+xn-1yn+xny1-x2y1-x3y2-⋯-xnyn-1-x1yn=12∑i=1n-1xiyi+1+xny1-∑i=1n-1xi+1yi-x1yn, (1)式中：S为锚位边界多边形的面积；(xi,yi)  (i=1,2,⋯,n )为顺时针排序的多边形顶点的经纬度坐标，其中xi为经度，yi为纬度．1.3节求取的锚位边界点列表是顺时针排序的，将其直接输入式(1)即可完成锚位面积的求取．依据锚位提取算法提取的锚位点是船舶AIS位置，并不是船舶的实际落锚位置．锚泊过程中，船舶围绕落锚点运动，运动范围受船舶自身长度、出链长度、水深等影响．本研究选取锚位边界多边形的形心作为所求船舶的落锚点的概位．求取边界的形心[17]的表达式为Cx=16S∑j=0n-1(xj+xj+1)(xjyj+1-xj+1yj)；(2)Cy=16S∑j=0n-1(yj+yj+1)(xjyj+1-xj+1yj)，(3)式中Cx和Cy分别为所求形心的经纬度坐标．2 实验分析2.1　实验数据为验证本文方法的可用性，本研究选取长江口锚地及其附近水域为实验区域．长江口锚地位于长江口定线制内，包括长江口NO.1~NO.3锚地，锚地形状规则如图7所示，图中：①为长江口NO.1锚地；②为长江口NO.2锚地；③为长江口NO.3锚地；品红色虚线圆为官方推荐锚位．长江口NO.1锚地位于长江口A警戒区东北侧，该锚地主要是供大型船舶和拟经长江口深水航道进口的船舶使用．长江口NO.2锚地位于长江口深水航道A警戒区东南侧，该锚地主要是供中小型船舶和拟经南槽航道进口的船舶使用．长江口NO.3锚地位于长江口A警戒区西北侧，该锚地主要是供拟经长江口深水航道进港船舶使用．采集该水域2019年4月1~2日的1 186艘船舶的3.168 95×105条AIS数据．原始船舶轨迹图如图8(a)所示．经去重、去噪、分段等预处理后得到3.161 00×105条有效AIS数据，预处理后的船舶轨迹图如图8(b)所示．10.13245/j.hust.239467.F007图7长江口锚地分布图10.13245/j.hust.239467.F008图82019年4月1~2日长江口锚地船舶轨迹图2.2　实验结果分析《海港总体设计规范》[18](以下简称《规范》)对船舶锚泊时所需的水域尺度给出了相关建议，指出船舶采用单锚系泊时，在理想状态下每个锚位所占水域为一圆形区域，其半径计算公式为R=L+3h+90，(4)式中：R为锚泊半径；L为设计船长；h为锚地水深．以式(4)为参照，对本文算法进行分析．首先，利用1.2节锚位识别算法提取船舶锚位，其中依据相关航海经验距离阈值设置为0.2 n mile，时间阈值设置为30 min，最终识别到锚泊船只200余艘次．然后，计算锚泊船舶的锚位面积与落锚点概位，并根据所求锚泊船舶数据形成锚泊船舶数据记录集，记录相关船舶的抛起锚时间、锚位面积、落锚点等信息，如表1所示，表中Sd为根据《规范》得到的锚位面积(锚地水深h=25 m)．10.13245/j.hust.239467.T001表12019年4月1~2日长江口锚地锚泊船舶信息表MMSI船长/m船宽/m吃水/m船舶类型开始时间2019-04-01结束时间2019-04-02落锚点概位S/(105m2)Sd/(105m2)3575730002293211.0散货00:01:4823:55:47122°38.111′E/31°9.328′N3.267 494.876 882320029692924510.5散货08:31:4004:04:35122°40.357′E/31°9.456′N3.274 206.561 1837458400090144.0油船00:03:0323:54:13122°41.906′E/31°1.939′N0.771 202.042 82413332130118185.5油船00:07:3410:01:45122°27.061′E/30°58.996′N1.102 432.516 07413526070164286.2集装箱15:10:2523:59:55122°27.624′E/31°11.823′N2.685 153.400 49413788000149218.5散杂货12:57:1703:34:39122°34.679′E/31°1.663′N2.093 333.097 48414313000190326.0散货00:01:4523:49:44122°42.815′E/31°2.845′N2.649 463.959 193705870003364613.6集装箱09:14:0723:58:56122°37.142′E/31°12.815′N2.153 857.885 42…………………………图9为2019年4月1~2日锚泊船舶的轨迹图，其中绿线代表船舶轨迹．以其中一艘巴拿马籍船舶SAKIZAYA DIAMOND轮(MMSI为357573000，船长为229 m，船宽为32 m，吃水为11.5 m，类型为散货船)为例进行详细说明．图10为该轮4月1~2日在长江口水域的历史轨迹，红色标签标注部分为船舶锚泊期间的轨迹．图11中绿色点为该船锚泊期间的轨迹点，蓝色多边形为本研究所提取的该船的锚位边界轮廓，红色点为本研究所求的该轮落锚点的概位．10.13245/j.hust.239467.F009图92019年4月1~2日长江口锚地锚泊船舶轨迹图10.13245/j.hust.239467.F010图10SAKIZAYA DIAMOND轮船舶轨迹10.13245/j.hust.239467.F011图11SAKIZAYA DIAMOND轮锚位面积计算表1统计结果显示，按照《规范》得到的锚位面积比根据本文算法计算而来的实际锚位面积要大．根据《规范》得到的锚位面积为船舶锚泊水域面积的推荐值，在此水域覆盖范围内，船舶是相对安全的．但是，占用面积过大也会降低锚地利用率．因此，当有锚地资源异常紧张的情况出现时，在气象水文条件允许，并留有足够的安全余量的情况下，可以考虑增加适当的锚位数量，降低船舶平均锚位面积占有量．在锚泊船舶的安全监控中，可以推荐值为安全阈值，当船舶实际锚位面积超过此值，则予以提醒存在走锚或碰撞的风险．表1统计结果还显示长江口锚地船舶平均锚泊时间较长，而这会大大降低锚地的运行效率，造成大量船舶积压在锚地，形成进不去、出不来的尴尬局面．在当前全球航运节能减排的大背景下，航运公司应合理制订船舶运营计划，合理使用经济航速，避免出现船舶锚地长时间等待的问题．3 结语本研究将数据挖掘算法与计算几何学相结合，提出一种锚位面积计算方法．选取长江口锚地水域的船舶AIS数据进行实例验证，结果表明：利用本文方法得到船舶锚位面积是可行的，该方法为锚地锚位的合理规划和锚泊船舶的安全监控提供数据支撑，对提高船舶锚地利用率和安全性具有积极意义．