王书献,孙永文,张胜茂,隋江华,朱文斌,杨胜龙,樊 伟
(1.大连海洋大学航海与船舶工程学院,辽宁大连 116023;2.浙江省海洋渔业资源可持续利用技术研究重点实验室,浙江省海洋水产研究所,浙江舟山 316021;3.农业农村部远洋与极地渔业创新重点实验室,中国水产科学研究院东海水产研究所,上海 200090)
船舶自动识别系统(automatic identification system,AIS)是由船载设备和岸基(星基)设施共同组成的一种数字系统。该系统将现代电子信息技术与网络技术结合,通过甚高频(very high frequency,VHF)频道向附近水域的岸台和船舶广播船位、船速及航向等船舶动态信息以及船名、吃水及危险货物运载状况等船舶静态资料[1-3]。国际海事组织发布的《国际海上人命安全公约》(International Convention for Safety of Life at Sea,SOLAS)要求航行于国际水域且总吨位在300以上的船舶,以及所有不论吨位大小的客船,均应安装AIS。AIS的广泛应用不仅对船舶间的信息沟通、碰撞避免以及内河航运的监管有着非常重要的作用,还为渔业发展带来了新的可能。
与传统陆基AIS相比,卫星AIS通过使用低轨小卫星或小卫星星座接收船载AIS站台发出的包含船舶静、动态数据的AIS信号,并将其转发到地面站进行分析及处理,从而实现更大范围甚至全球海域的舰船监视。目前,国内外学者对AIS数据进行了分析挖掘,已在渔场分析、非法捕捞行为监控、海上交通安全等方面取得了一系列成果[4-6]。周海等[7]提出了一种从AIS数据中提取出船舶运动模式的方法,丰富了船舶轨迹聚类的模式特征;COURIE IA等[8]提出利用AIS数据设计微型卫星成像系统,预防印度尼西亚非法捕捞;GREIG N C等[9]从海上交通安全与生物安全的角度,提出使用卫星AIS数据分析美国华盛顿州船舶航行速度,并将其作为鲸鱼类生物与船舶碰撞分析的依据。但是,卫星AIS数据在船舶时空性分析(尤其是渔船位置分析与可视化)方面依然存在很大的潜在价值等待挖掘。
国家卫星海洋应用中心接收有我国卫星AIS数据。该系统每天接收的AIS数据约有50万条,分别单独存储于以日期命名的文件夹中。本文利用的卫星AIS数据字段包含船舶唯一标识符(maritimemobile service identify,MMSI)、经度、纬度、航速和航向等,并以热力图的形式展示船舶空间分布情况。
热力图是一个以颜色变化来显示数据的矩阵。具体来说,通过密度函数对输入数据进行可视化,并以特殊高亮的形式在底面图层显示。一般而言,从数据类型看,热力图是用2个连续字段确定数值点的位置;从功能看,热力图主要用于展示数据的分布情况。热力图凭借其直观性和便捷性,逐渐在交通、环境等各个领域中得到了广泛的应用。詹显军等[10]以热力图展示地铁车站在施工过程中整体变形的实时分布情况,提出了“基于热力图的地铁车站变形可视化分析”;吕玉嫦等[11]提出热力图在气象网站建设中的应用,美观、直接地展示天气情况。但是热力图制图技术在根据卫星AIS数据监测渔船位置方面的应用相对比较缺乏。以热力图的形式展示海量的渔船位置信息,可以将抽象的经纬数据可视化在一张地图中,便于观察及进一步分析、决策[12-15]。生成热力图的方法有多种,其中,使用arcpy站点包的方法不仅可以实现生成过程的自动化,还方便后期修改和完善热力图。
arcpy是ArcGIS的一个站点包。该站点包以arcgisscripting模块为基础。它的目的是通过Python语言以高效实用的方式执行数据的管理、制作地图的自动化以及地理数据的分析。arcpy站点包的接口和函数,均符合习惯性的思维模式。其设计都以精简实用为主导,避免了像其他高级编程语言那样严格的定义语法、变量类型等格式化的操作,符合Python的语法特点[16-20]。
本文基于卫星AIS数据,利用Python程序对数据进行筛选、分类,并提取出有效信息。最后将有效信息导入ArcGIS地图工程,利用arcpy程序包,按照日期分别自动渲染成热力图,最终形成专题热力图。利用该专题图,监测和统计渔船在各海域的分布情况,了解捕捞渔船动态变化,分析渔场时空特征,辅助渔业捕捞决策。
AIS数据存储于国家卫星海洋应用中心FTP(ftp.nsoas.org.cn)。本文使用的卫星AIS数据均从该网站下载。每天的AIS数据存储在一个文件夹内,包含后缀分别为jpg、l1b、xml的3个文件。命名规则为“MUL_OPEN_AIS_L1B_”拼接上“日期”及“_SSS”,例如,2020年5月14日的数据文件夹中分别包含“MUL_OPEN_AIS_L1B_20200514_SSS.jpg”“MUL_OPEN_AIS_L1B_20200514_SSS.l1b”和“MUL_OPEN_AIS_L1B_20200514_SSS.meta.xml”3个文件。其中jpg文件为缩略图文件,l1b文件是以CSV形式存储的AIS数据,xml为元数据文件。
根据海洋观测卫星地面系统AIS 0级和1级产品数据格式(海洋一号C/D星和海洋二号B/C星搭载的AIS载荷0级和1级产品)中的消息类型,l1b文件中按照不同格式存储有船舶唯一标识符MMSI、消息类型、经度、纬度、航速、航向、呼号和名称等信息。该文件将这些数据以CSV形式存储,每条消息报文占一行,每条消息报文的字段值之间以半角逗号分隔。
由于AIS数据每天更新,手动下载虽然可以完成目标,但是十分费时费力。利用Python程序,模拟手动获取数据的过程,可提高工作效率。定义一个MyFtp类,利用ftplib库,编写初始化方法、登录方法、上传文件方法和下载文件方法等。这些方法编写完成后,一个简单通用的FTP类即完成。在出现具体的FTP需求时,可以编写新类继承自该类(图1)。在本研究中,编写AISftp类继承自MyFTP类,用于实现卫星AIS数据的读取和下载。
利用自定义的FTP相关类,配置好FTP信息并完成FTP下载(图2)。首先,配置好主机、端口号、用户名和密码等信息,登录FTP;其次,利用look_for_files()方法在FTP指定目录中找到所有后缀为“.l1b”的文件,最后将其依次下载到本地磁盘的指定目录中。
本研究中测试FTP客户端连接使用FlashFXP 5 软件,开发环境的操作系统为Windows 10专业版,开发语言为Python 3.8,Python开发IDE为PyCharm。热力图制作使用ArcGIS,由于arcpy包不支持由Python Package Index(PYPI)直接安装,且ArcGIS 10.7推荐使用内置的Python 2,因此本文中Arcpy的运行环境为ArcGIS内置的Python 2.7。
图1 FTP相关类的结构Fig.1 Structure of FTP related classes
图2 FTP下载流程Fig.2 FTP download process
卫星AIS数据有27种标准报文,其报文字段并不完全一致。对本研究而言,有价值的报文应包含有日期、时间、船舶唯一标识符、经度和纬度等字段。在该27类报文中,符合上述条件的报文类型有:MESSAGE 1、MESSAGE 2、MESSAGE 3、MESSAGE 4、MESSAGE 11、MESSAGE 18、MESSAGE 19和MESSAGE 27。其中MESSAGE 1和MESSAGE 2的存储格式完全一致,MESSAGE 4和MESSAGE 11的存储格式完全一致。
基于以上分析,可以遍历从FTP下载的l1b数据,按照消息类型分别进行处理。若消息类型为上述分析中含有关键信息的类型,则存入对应的数据库中;否则放弃该消息报文,直接处理下一条。
使用MySQL数据库,创建一个AIS库,按照不同消息类型的存储格式分别创建好对应的数据表。以MESSAGE 1和MESSAGE 2为例,从“十二五”时期海洋观测卫星地面系统AIS 0和AIS 1级产品数据格式文档中查找到MESSAGE 1和MESSAGE 2的报文格式,如表1所示。
按照该格式在数据库中新建一个MESSAGE 1and MESSAGE 2表,每个报文字段作为一个数据库字段。其他消息类型可仿照MESSAGE 1和MESSAGE 2的操作方法分别创建数据库表。待下载的卫星AIS数据遍历完成后,所有有效消息报文即按照消息类型分别存储在数据库的不同表中。
MMSI是船舶无线电通信系统在其无线电信道上发送的、能独特识别各类台站和成组呼叫台站的一列9位数字码[21-22]。由于该数字码可以唯一标识一艘船舶,故可以收集所有渔船的MMSI,对比报文中MMSI与收集到的渔船MMSI,以区分渔船及非渔船。若报文中的MMSI存在于收集到的渔船MMSI列表中,则该条卫星AIS数据为渔船数据,否则为非渔船数据。
从中国远洋VMS系统等平台收集国内外渔船的MMSI记录共25 995条(中国远洋VMS系统里包括船舶MMSI、类型和吨位信息。历史共计3 260条,没有MMSI号记录20条,NULL记录7条。将该数据补充到global fishing watch里面,共计25 995条渔船MMSI)。将该25 995条MMSI写入fishing_boat.csv文件中,便于程序调用。
简单统计处理后的数据,得到一个明显的基础性结论:每日观测到的渔船数量远低于非渔船的数量。饼状图可以直观地展示其比例关系。以5月3日和5月4日为例:5月3日监测到的渔船数量为2 975艘,占比12.66%,非渔船数量20 520艘,占比87.34%(图3-a);5月4日监测到的渔船数量为2 971艘,占比12.55%,非渔船数量20 701艘,占比87.45%(图3-b)。
表1 MESSAGE 1和MESSAGE 2消息格式Tab.1 M essage format of MESSAGE 1 and MESSAGE 2
图3 部分船舶类型比例Fig.3 Proportion of some ship types
使用ArcGIS软件,新建一个地图文件,选择WorldRobinson模板作为热力图模板的底面图层。再将前文中获取的数据文件导入。设置数据坐标系为“WGS-1984”(与卫星AIS报文数据坐标系一致)。将每个渔船位置(经纬度作为输入)作为一个点,在地图上方展示出来(图4)。所得散点图可以清晰展示渔船分布状态。
图4 5月3日渔船散点图Fig.4 Scatter plot of fishing boat on M ay 3rd
点密度分析工具使用输入点要素或线要素来计算感兴趣区域内的密度地图。该工具用于计算每个输出栅格像元周围的点要素密度[18]。在每个栅格像元中心的周围定义一个邻域,该邻域的半径大小可以设定:半径参数值越大,生成的密度栅格概化程度便越高;值越小,生成的栅格所显示信息越详细。将该邻域内点的数量相加,然后除以邻域面积,即得到点要素的密度[19],即点要素(x,y)的密度计算公式如下。
式(1)中,(x,y)表示横坐标为x,纵坐标为y的点要素,∂(x,y)表示点要素(x,y)的某个邻域,S∂(x,y)表示点要素(x,y)的邻域面积,ρ∂(x,y)表示点要素(x,y)的密度值。在一定范围内扩大或缩小半径值,并不会明显改变密度值。因为密度值为邻域内点数与邻域面积的比值,而邻域内点数量的变化是伴随着邻域面积变化的。因此,邻域面积的大小实际影响的是热力图的概化程度。点密度分析工具根据式(1)分析所有数据点,并以不同颜色表示不同密度值高亮显示(图5)。为达到较好的成图效果,本文点密度分析中邻域半径的确定方式为:首先,计算出输出空间参考中长度、宽度的最小值作为基数;再将该基数乘以系数k(本文k取1/30)作为邻域半径。选取输出空间参考中长度、宽度的最小值作为邻域半径确定的基数是因为当输出空间参考长度、宽度相差较大时,若选取最大值作为邻域半径确定基数,会导致热力图在最小值方向上概化程度过高,无法精细反映热力分布状态。本文中系数k取ArcGIS官方推荐的1/30,成图结果表明,该推荐系数可以较好地反映渔船分布状态。
将渔船位置信息作为点要素,制作热力图模板,可以在直观观测渔船分布的基础上对渔船未来分布趋势作出预测。
图5 点密度分析工作流Fig.5 Point density analysis work-flow
arcpy是一个Python站点包,可提供一系列实用高效的方法,执行地理数据分析、数据转换、数据管理和地图自动化。在pycharm 中编写Python程序,并将其运行环境配置为ArcGIS内置的Python 2.7环境,即可直接在程序中导入arcpy程序包并使用。利用arcpy中定义的类及其方法,可以在不打开ArcGIS软件的情况下,完成对地图文件的编辑和导出。使用arcpy程序包,为每一天的数据生成热力图,最终形成专题热力图。
首先,编写生成某一天的散点图、热力图和含散点热力图的方法(图6-a),将日期作为该方法的参数。
每天的热力图生成方法完成后,在主程序中遍历获取到的卫星AIS数据日期,依次传入该方法(图6-b)。
程序执行完毕后,每个日期分别生成3张图。以5月3日为例,生成:散点图、热力图(图7-a)和含散点热力图(图7-b)。
本文使用国家卫星海洋应用中心FTP(ftp.nsoas.org.cn)存储的卫星AIS数据,利用Python程序及相关类库,实现了数据下载、数据分类、数据分析和专题数据自动化制图等,并以热力图直观展示渔船位置信息,为渔业部门的数据统计、未来决策等提供数据参考,推动渔业发展。
从FTP下载、数据解析、数据筛选到数据绘制,Python程序取代人工完成了大量的重复工作。相较于人工操作,本课题在实现过程中节约了大量的时间。但是,可以改进的空间仍然很大。例如,实验过程中数据下载部分耗时较久,这是因为FTP端设置了速度限制。可以在脚本中加入多线程,以多个连接同时下载,理论上可以大幅提升下载速度;在数据分类分析中,本文使用了MySQL数据库,频繁的数据库插入操作十分耗时,设计本地硬盘处理方案可以加快程序整体处理速度。
图6 程序流程图Fig.6 Program flow
点密度分析是本文热力图制作的基础,热力图不同颜色表示不同的密度值。密度在物理学上的定义为:单位体积内物体的质量。本文中密度值的物理意义表现为单位面积内渔船的数量。在点密度分析研究中,空间尺度(邻域)的确定往往会影响到热力图是否美观、直观与合理[22-23]。邻域半径的设置会影响到热力图的概化程度。例如,图8-a中存在一组测试数据,分别以半径为8.278 cm、2.54 cm和25.4 cm的邻域对该组数据进行点密度分析,形成了图8-b、图8-c和图8-d。结果表明,若邻域半径设置过小,热力图概化程度过低,易出现类似于二值化的热力图。图8-c中点几乎与图8-a完全一致,虽然能够精细地表示每个数据的位置,但是放弃了热力图直观表示数据分布情况的目标;若邻域半径设置过大,则热力图概化程度过高,某些区域的船舶密度表达不清晰。为了避免上述2种极端情况,本文采用的邻域半径确定方法为:找出输出数据区域长度、宽度中的最小值,将该值除以30作为点密度分析邻域的半径。实验结果表明,该方案可以达到预期目标。
本文利用Python程序对AIS卫星监测的船舶信息进行筛选、分类和可视化,最终形成了专题热力图。将复杂的渔船分布信息以热力图的形式在地图中显示,不仅节约了人工分析AIS数据的时间,还让抽象数据变得更直观。研究表明:
(1)AIS数据的下载、分类、分析、专题图自动制图等流程可以由Python程序完成。节约大量人工成本。
(2)点密度分析中,空间尺度的设定会通过数据概化程度影响热力图制图效果。若空间尺度设置过小,则热力图会出现二值化现象,与原始数据散点一致,无热力参考价值;若空间尺度设置过大,则热力图概化程度过高,无法精确体现热力分布。以输出区域长度、宽度中最小值的1/30作为空间尺度半径,可以避免上述极端情况,是一个较好的折衷方案。
图7 程序输出热力图Fig.7 Heatm aps by program
图8 邻域半径的确定Fig.8 Determ ination of neighborhood radius
(3)从专题热力图结果看,每天的渔船船位分布并不一致,但热力分布规律仍有规律可循。例如:与阿根廷接壤的麦哲伦海峡、南非伊丽莎白港附近海域、与澳大利亚接壤的塔斯曼海域、新西兰附近海域、与俄罗斯东部接壤的鄂霍次克海域等区域具有较大的分布密度。尤其是与阿根廷接壤的麦哲伦海峡,在各个日期内的热力值均处于最高状态。
(4)利用程序生成专题热力图,可从热力图中发现渔船分布规律,为渔业企业及政府部门未来规划、决策提供直观数据参考。