王書獻(xiàn)，孫永文，張勝茂，隋江華，朱文斌，楊勝龍，樊 偉

（1.大連海洋大學(xué)航海與船舶工程學(xué)院，遼寧大連 116023；2.浙江省海洋漁業(yè)資源可持續(xù)利用技術(shù)研究重點實驗室，浙江省海洋水產(chǎn)研究所，浙江舟山 316021；3.農(nóng)業(yè)農(nóng)村部遠(yuǎn)洋與極地漁業(yè)創(chuàng)新重點實驗室，中國水產(chǎn)科學(xué)研究院東海水產(chǎn)研究所，上海 200090）

船舶自動識別系統(tǒng)（automatic identification system，AIS）是由船載設(shè)備和岸基（星基）設(shè)施共同組成的一種數(shù)字系統(tǒng)。該系統(tǒng)將現(xiàn)代電子信息技術(shù)與網(wǎng)絡(luò)技術(shù)結(jié)合，通過甚高頻（very high frequency，VHF）頻道向附近水域的岸臺和船舶廣播船位、船速及航向等船舶動態(tài)信息以及船名、吃水及危險貨物運載狀況等船舶靜態(tài)資料［1－3］。國際海事組織發(fā)布的《國際海上人命安全公約》（International Convention for Safety of Life at Sea，SOLAS）要求航行于國際水域且總噸位在300以上的船舶，以及所有不論噸位大小的客船，均應(yīng)安裝AIS。AIS的廣泛應(yīng)用不僅對船舶間的信息溝通、碰撞避免以及內(nèi)河航運的監(jiān)管有著非常重要的作用，還為漁業(yè)發(fā)展帶來了新的可能。

與傳統(tǒng)陸基AIS相比，衛(wèi)星AIS通過使用低軌小衛(wèi)星或小衛(wèi)星星座接收船載AIS站臺發(fā)出的包含船舶靜、動態(tài)數(shù)據(jù)的AIS信號，并將其轉(zhuǎn)發(fā)到地面站進(jìn)行分析及處理，從而實現(xiàn)更大范圍甚至全球海域的艦船監(jiān)視。目前，國內(nèi)外學(xué)者對AIS數(shù)據(jù)進(jìn)行了分析挖掘，已在漁場分析、非法捕撈行為監(jiān)控、海上交通安全等方面取得了一系列成果［4－6］。周海等［7］提出了一種從AIS數(shù)據(jù)中提取出船舶運動模式的方法，豐富了船舶軌跡聚類的模式特征；COURIE IA等［8］提出利用AIS數(shù)據(jù)設(shè)計微型衛(wèi)星成像系統(tǒng)，預(yù)防印度尼西亞非法捕撈；GREIG N C等［9］從海上交通安全與生物安全的角度，提出使用衛(wèi)星AIS數(shù)據(jù)分析美國華盛頓州船舶航行速度，并將其作為鯨魚類生物與船舶碰撞分析的依據(jù)。但是，衛(wèi)星AIS數(shù)據(jù)在船舶時空性分析（尤其是漁船位置分析與可視化）方面依然存在很大的潛在價值等待挖掘。

國家衛(wèi)星海洋應(yīng)用中心接收有我國衛(wèi)星AIS數(shù)據(jù)。該系統(tǒng)每天接收的AIS數(shù)據(jù)約有50萬條，分別單獨存儲于以日期命名的文件夾中。本文利用的衛(wèi)星AIS數(shù)據(jù)字段包含船舶唯一標(biāo)識符（maritimemobile service identify，MMSI）、經(jīng)度、緯度、航速和航向等，并以熱力圖的形式展示船舶空間分布情況。

熱力圖是一個以顏色變化來顯示數(shù)據(jù)的矩陣。具體來說，通過密度函數(shù)對輸入數(shù)據(jù)進(jìn)行可視化，并以特殊高亮的形式在底面圖層顯示。一般而言，從數(shù)據(jù)類型看，熱力圖是用2個連續(xù)字段確定數(shù)值點的位置；從功能看，熱力圖主要用于展示數(shù)據(jù)的分布情況。熱力圖憑借其直觀性和便捷性，逐漸在交通、環(huán)境等各個領(lǐng)域中得到了廣泛的應(yīng)用。詹顯軍等［10］以熱力圖展示地鐵車站在施工過程中整體變形的實時分布情況，提出了“基于熱力圖的地鐵車站變形可視化分析”；呂玉嫦等［11］提出熱力圖在氣象網(wǎng)站建設(shè)中的應(yīng)用，美觀、直接地展示天氣情況。但是熱力圖制圖技術(shù)在根據(jù)衛(wèi)星AIS數(shù)據(jù)監(jiān)測漁船位置方面的應(yīng)用相對比較缺乏。以熱力圖的形式展示海量的漁船位置信息，可以將抽象的經(jīng)緯數(shù)據(jù)可視化在一張地圖中，便于觀察及進(jìn)一步分析、決策［12－15］。生成熱力圖的方法有多種，其中，使用arcpy站點包的方法不僅可以實現(xiàn)生成過程的自動化，還方便后期修改和完善熱力圖。

arcpy是ArcGIS的一個站點包。該站點包以arcgisscripting模塊為基礎(chǔ)。它的目的是通過Python語言以高效實用的方式執(zhí)行數(shù)據(jù)的管理、制作地圖的自動化以及地理數(shù)據(jù)的分析。arcpy站點包的接口和函數(shù)，均符合習(xí)慣性的思維模式。其設(shè)計都以精簡實用為主導(dǎo)，避免了像其他高級編程語言那樣嚴(yán)格的定義語法、變量類型等格式化的操作，符合Python的語法特點［16－20］。

本文基于衛(wèi)星AIS數(shù)據(jù)，利用Python程序?qū)?shù)據(jù)進(jìn)行篩選、分類，并提取出有效信息。最后將有效信息導(dǎo)入ArcGIS地圖工程，利用arcpy程序包，按照日期分別自動渲染成熱力圖，最終形成專題熱力圖。利用該專題圖，監(jiān)測和統(tǒng)計漁船在各海域的分布情況，了解捕撈漁船動態(tài)變化，分析漁場時空特征，輔助漁業(yè)捕撈決策。

1 數(shù)據(jù)統(tǒng)計及處理

1.1 數(shù)據(jù)獲取

AIS數(shù)據(jù)存儲于國家衛(wèi)星海洋應(yīng)用中心FTP（ftp.nsoas.org.cn）。本文使用的衛(wèi)星AIS數(shù)據(jù)均從該網(wǎng)站下載。每天的AIS數(shù)據(jù)存儲在一個文件夾內(nèi)，包含后綴分別為jpg、l1b、xml的3個文件。命名規(guī)則為“MUL＿OPEN＿AIS＿L1B＿”拼接上“日期”及“＿SSS”，例如，2020年5月14日的數(shù)據(jù)文件夾中分別包含“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數(shù)據(jù)，xml為元數(shù)據(jù)文件。

根據(jù)海洋觀測衛(wèi)星地面系統(tǒng)AIS 0級和1級產(chǎn)品數(shù)據(jù)格式（海洋一號C／D星和海洋二號B／C星搭載的AIS載荷0級和1級產(chǎn)品）中的消息類型，l1b文件中按照不同格式存儲有船舶唯一標(biāo)識符MMSI、消息類型、經(jīng)度、緯度、航速、航向、呼號和名稱等信息。該文件將這些數(shù)據(jù)以CSV形式存儲，每條消息報文占一行，每條消息報文的字段值之間以半角逗號分隔。

由于AIS數(shù)據(jù)每天更新，手動下載雖然可以完成目標(biāo)，但是十分費時費力。利用Python程序，模擬手動獲取數(shù)據(jù)的過程，可提高工作效率。定義一個MyFtp類，利用ftplib庫，編寫初始化方法、登錄方法、上傳文件方法和下載文件方法等。這些方法編寫完成后，一個簡單通用的FTP類即完成。在出現(xiàn)具體的FTP需求時，可以編寫新類繼承自該類（圖1）。在本研究中，編寫AISftp類繼承自MyFTP類，用于實現(xiàn)衛(wèi)星AIS數(shù)據(jù)的讀取和下載。

利用自定義的FTP相關(guān)類，配置好FTP信息并完成FTP下載（圖2）。首先，配置好主機(jī)、端口號、用戶名和密碼等信息，登錄FTP；其次，利用look＿for＿files（）方法在FTP指定目錄中找到所有后綴為“.l1b”的文件，最后將其依次下載到本地磁盤的指定目錄中。

本研究中測試FTP客戶端連接使用FlashFXP 5 軟件，開發(fā)環(huán)境的操作系統(tǒng)為Windows 10專業(yè)版，開發(fā)語言為Python 3.8，Python開發(fā)IDE為PyCharm。熱力圖制作使用ArcGIS，由于arcpy包不支持由Python Package Index（PYPI）直接安裝，且ArcGIS 10.7推薦使用內(nèi)置的Python 2，因此本文中Arcpy的運行環(huán)境為ArcGIS內(nèi)置的Python 2.7。

圖1 FTP相關(guān)類的結(jié)構(gòu)Fig.1 Structure of FTP related classes

圖2 FTP下載流程Fig.2 FTP download process

1.2 數(shù)據(jù)處理

衛(wèi)星AIS數(shù)據(jù)有27種標(biāo)準(zhǔn)報文，其報文字段并不完全一致。對本研究而言，有價值的報文應(yīng)包含有日期、時間、船舶唯一標(biāo)識符、經(jīng)度和緯度等字段。在該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數(shù)據(jù)，按照消息類型分別進(jìn)行處理。若消息類型為上述分析中含有關(guān)鍵信息的類型，則存入對應(yīng)的數(shù)據(jù)庫中；否則放棄該消息報文，直接處理下一條。

使用MySQL數(shù)據(jù)庫，創(chuàng)建一個AIS庫，按照不同消息類型的存儲格式分別創(chuàng)建好對應(yīng)的數(shù)據(jù)表。以MESSAGE 1和MESSAGE 2為例，從“十二五”時期海洋觀測衛(wèi)星地面系統(tǒng)AIS 0和AIS 1級產(chǎn)品數(shù)據(jù)格式文檔中查找到MESSAGE 1和MESSAGE 2的報文格式，如表1所示。

按照該格式在數(shù)據(jù)庫中新建一個MESSAGE 1and MESSAGE 2表，每個報文字段作為一個數(shù)據(jù)庫字段。其他消息類型可仿照MESSAGE 1和MESSAGE 2的操作方法分別創(chuàng)建數(shù)據(jù)庫表。待下載的衛(wèi)星AIS數(shù)據(jù)遍歷完成后，所有有效消息報文即按照消息類型分別存儲在數(shù)據(jù)庫的不同表中。

MMSI是船舶無線電通信系統(tǒng)在其無線電信道上發(fā)送的、能獨特識別各類臺站和成組呼叫臺站的一列9位數(shù)字碼［21－22］。由于該數(shù)字碼可以唯一標(biāo)識一艘船舶，故可以收集所有漁船的MMSI，對比報文中MMSI與收集到的漁船MMSI，以區(qū)分漁船及非漁船。若報文中的MMSI存在于收集到的漁船MMSI列表中，則該條衛(wèi)星AIS數(shù)據(jù)為漁船數(shù)據(jù)，否則為非漁船數(shù)據(jù)。

從中國遠(yuǎn)洋VMS系統(tǒng)等平臺收集國內(nèi)外漁船的MMSI記錄共25 995條（中國遠(yuǎn)洋VMS系統(tǒng)里包括船舶MMSI、類型和噸位信息。歷史共計3 260條，沒有MMSI號記錄20條，NULL記錄7條。將該數(shù)據(jù)補(bǔ)充到global fishing watch里面，共計25 995條漁船MMSI）。將該25 995條MMSI寫入fishing＿boat.csv文件中，便于程序調(diào)用。

1.3 數(shù)據(jù)統(tǒng)計

簡單統(tǒng)計處理后的數(shù)據(jù)，得到一個明顯的基礎(chǔ)性結(jié)論：每日觀測到的漁船數(shù)量遠(yuǎn)低于非漁船的數(shù)量。餅狀圖可以直觀地展示其比例關(guān)系。以5月3日和5月4日為例：5月3日監(jiān)測到的漁船數(shù)量為2 975艘，占比12.66%，非漁船數(shù)量20 520艘，占比87.34%（圖3－a）；5月4日監(jiān)測到的漁船數(shù)量為2 971艘，占比12.55%，非漁船數(shù)量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

2 熱力圖生成

2.1 熱力圖模板制作

使用ArcGIS軟件，新建一個地圖文件，選擇WorldRobinson模板作為熱力圖模板的底面圖層。再將前文中獲取的數(shù)據(jù)文件導(dǎo)入。設(shè)置數(shù)據(jù)坐標(biāo)系為“WGS-1984”（與衛(wèi)星AIS報文數(shù)據(jù)坐標(biāo)系一致）。將每個漁船位置（經(jīng)緯度作為輸入）作為一個點，在地圖上方展示出來（圖4）。所得散點圖可以清晰展示漁船分布狀態(tài)。

圖4 5月3日漁船散點圖Fig.4 Scatter plot of fishing boat on M ay 3rd

點密度分析工具使用輸入點要素或線要素來計算感興趣區(qū)域內(nèi)的密度地圖。該工具用于計算每個輸出柵格像元周圍的點要素密度［18］。在每個柵格像元中心的周圍定義一個鄰域，該鄰域的半徑大小可以設(shè)定：半徑參數(shù)值越大，生成的密度柵格概化程度便越高；值越小，生成的柵格所顯示信息越詳細(xì)。將該鄰域內(nèi)點的數(shù)量相加，然后除以鄰域面積，即得到點要素的密度［19］，即點要素（x，y）的密度計算公式如下。

式（1）中，（x，y）表示橫坐標(biāo)為x，縱坐標(biāo)為y的點要素，?（x，y）表示點要素（x，y）的某個鄰域，S?（x，y）表示點要素（x，y）的鄰域面積，ρ?（x，y）表示點要素（x，y）的密度值。在一定范圍內(nèi)擴(kuò)大或縮小半徑值，并不會明顯改變密度值。因為密度值為鄰域內(nèi)點數(shù)與鄰域面積的比值，而鄰域內(nèi)點數(shù)量的變化是伴隨著鄰域面積變化的。因此，鄰域面積的大小實際影響的是熱力圖的概化程度。點密度分析工具根據(jù)式（1）分析所有數(shù)據(jù)點，并以不同顏色表示不同密度值高亮顯示（圖5）。為達(dá)到較好的成圖效果，本文點密度分析中鄰域半徑的確定方式為：首先，計算出輸出空間參考中長度、寬度的最小值作為基數(shù)；再將該基數(shù)乘以系數(shù)k（本文k取1／30）作為鄰域半徑。選取輸出空間參考中長度、寬度的最小值作為鄰域半徑確定的基數(shù)是因為當(dāng)輸出空間參考長度、寬度相差較大時，若選取最大值作為鄰域半徑確定基數(shù)，會導(dǎo)致熱力圖在最小值方向上概化程度過高，無法精細(xì)反映熱力分布狀態(tài)。本文中系數(shù)k取ArcGIS官方推薦的1／30，成圖結(jié)果表明，該推薦系數(shù)可以較好地反映漁船分布狀態(tài)。

將漁船位置信息作為點要素，制作熱力圖模板，可以在直觀觀測漁船分布的基礎(chǔ)上對漁船未來分布趨勢作出預(yù)測。

圖5 點密度分析工作流Fig.5 Point density analysis work-flow

2.2 自動生成

arcpy是一個Python站點包，可提供一系列實用高效的方法，執(zhí)行地理數(shù)據(jù)分析、數(shù)據(jù)轉(zhuǎn)換、數(shù)據(jù)管理和地圖自動化。在pycharm 中編寫Python程序，并將其運行環(huán)境配置為ArcGIS內(nèi)置的Python 2.7環(huán)境，即可直接在程序中導(dǎo)入arcpy程序包并使用。利用arcpy中定義的類及其方法，可以在不打開ArcGIS軟件的情況下，完成對地圖文件的編輯和導(dǎo)出。使用arcpy程序包，為每一天的數(shù)據(jù)生成熱力圖，最終形成專題熱力圖。

首先，編寫生成某一天的散點圖、熱力圖和含散點熱力圖的方法（圖6－a），將日期作為該方法的參數(shù)。

每天的熱力圖生成方法完成后，在主程序中遍歷獲取到的衛(wèi)星AIS數(shù)據(jù)日期，依次傳入該方法（圖6－b）。

程序執(zhí)行完畢后，每個日期分別生成3張圖。以5月3日為例，生成：散點圖、熱力圖（圖7－a）和含散點熱力圖（圖7－b）。

3 討論及結(jié)論

本文使用國家衛(wèi)星海洋應(yīng)用中心FTP（ftp.nsoas.org.cn）存儲的衛(wèi)星AIS數(shù)據(jù)，利用Python程序及相關(guān)類庫，實現(xiàn)了數(shù)據(jù)下載、數(shù)據(jù)分類、數(shù)據(jù)分析和專題數(shù)據(jù)自動化制圖等，并以熱力圖直觀展示漁船位置信息，為漁業(yè)部門的數(shù)據(jù)統(tǒng)計、未來決策等提供數(shù)據(jù)參考，推動漁業(yè)發(fā)展。

從FTP下載、數(shù)據(jù)解析、數(shù)據(jù)篩選到數(shù)據(jù)繪制，Python程序取代人工完成了大量的重復(fù)工作。相較于人工操作，本課題在實現(xiàn)過程中節(jié)約了大量的時間。但是，可以改進(jìn)的空間仍然很大。例如，實驗過程中數(shù)據(jù)下載部分耗時較久，這是因為FTP端設(shè)置了速度限制。可以在腳本中加入多線程，以多個連接同時下載，理論上可以大幅提升下載速度；在數(shù)據(jù)分類分析中，本文使用了MySQL數(shù)據(jù)庫，頻繁的數(shù)據(jù)庫插入操作十分耗時，設(shè)計本地硬盤處理方案可以加快程序整體處理速度。

圖6 程序流程圖Fig.6 Program flow

點密度分析是本文熱力圖制作的基礎(chǔ)，熱力圖不同顏色表示不同的密度值。密度在物理學(xué)上的定義為：單位體積內(nèi)物體的質(zhì)量。本文中密度值的物理意義表現(xiàn)為單位面積內(nèi)漁船的數(shù)量。在點密度分析研究中，空間尺度（鄰域）的確定往往會影響到熱力圖是否美觀、直觀與合理［22－23］。鄰域半徑的設(shè)置會影響到熱力圖的概化程度。例如，圖8－a中存在一組測試數(shù)據(jù)，分別以半徑為8.278 cm、2.54 cm和25.4 cm的鄰域?qū)υ摻M數(shù)據(jù)進(jìn)行點密度分析，形成了圖8－b、圖8－c和圖8－d。結(jié)果表明，若鄰域半徑設(shè)置過小，熱力圖概化程度過低，易出現(xiàn)類似于二值化的熱力圖。圖8－c中點幾乎與圖8－a完全一致，雖然能夠精細(xì)地表示每個數(shù)據(jù)的位置，但是放棄了熱力圖直觀表示數(shù)據(jù)分布情況的目標(biāo)；若鄰域半徑設(shè)置過大，則熱力圖概化程度過高，某些區(qū)域的船舶密度表達(dá)不清晰。為了避免上述2種極端情況，本文采用的鄰域半徑確定方法為：找出輸出數(shù)據(jù)區(qū)域長度、寬度中的最小值，將該值除以30作為點密度分析鄰域的半徑。實驗結(jié)果表明，該方案可以達(dá)到預(yù)期目標(biāo)。

本文利用Python程序?qū)IS衛(wèi)星監(jiān)測的船舶信息進(jìn)行篩選、分類和可視化，最終形成了專題熱力圖。將復(fù)雜的漁船分布信息以熱力圖的形式在地圖中顯示，不僅節(jié)約了人工分析AIS數(shù)據(jù)的時間，還讓抽象數(shù)據(jù)變得更直觀。研究表明：

（1）AIS數(shù)據(jù)的下載、分類、分析、專題圖自動制圖等流程可以由Python程序完成。節(jié)約大量人工成本。

（2）點密度分析中，空間尺度的設(shè)定會通過數(shù)據(jù)概化程度影響熱力圖制圖效果。若空間尺度設(shè)置過小，則熱力圖會出現(xiàn)二值化現(xiàn)象，與原始數(shù)據(jù)散點一致，無熱力參考價值；若空間尺度設(shè)置過大，則熱力圖概化程度過高，無法精確體現(xiàn)熱力分布。以輸出區(qū)域長度、寬度中最小值的1／30作為空間尺度半徑，可以避免上述極端情況，是一個較好的折衷方案。

圖7 程序輸出熱力圖Fig.7 Heatm aps by program

圖8 鄰域半徑的確定Fig.8 Determ ination of neighborhood radius

（3）從專題熱力圖結(jié)果看，每天的漁船船位分布并不一致，但熱力分布規(guī)律仍有規(guī)律可循。例如：與阿根廷接壤的麥哲倫海峽、南非伊麗莎白港附近海域、與澳大利亞接壤的塔斯曼海域、新西蘭附近海域、與俄羅斯東部接壤的鄂霍次克海域等區(qū)域具有較大的分布密度。尤其是與阿根廷接壤的麥哲倫海峽，在各個日期內(nèi)的熱力值均處于最高狀態(tài)。

（4）利用程序生成專題熱力圖，可從熱力圖中發(fā)現(xiàn)漁船分布規(guī)律，為漁業(yè)企業(yè)及政府部門未來規(guī)劃、決策提供直觀數(shù)據(jù)參考。