楊家軒， 吳京霖， 姜大鵬

(大連海事大學 a.航海學院; b.遼寧省航海安全保障重點實驗室， 遼寧 大連 116026)

隨著計算機技術(shù)和通信技術(shù)的快速發(fā)展，各領(lǐng)域的信息暴增，人類社會已進入大數(shù)據(jù)高速發(fā)展的時代，如何從大量的數(shù)據(jù)中發(fā)現(xiàn)潛在的特征，把握區(qū)域事件的發(fā)展規(guī)律已成為相關(guān)學者研究的焦點。近年來,航運信息的發(fā)展為海上交通領(lǐng)域積累了海量的數(shù)據(jù)，其中以船舶自動識別系統(tǒng)(Automatic Identification System, AIS)為主的船舶歷史數(shù)據(jù)扮演著重要角色。AIS的廣泛應用不僅提高船舶航行的安全性,而且其數(shù)據(jù)能較好地反映港口水域船舶海上交通信息，分析大量的位置信息、航線數(shù)據(jù)和AIS 數(shù)據(jù)等，尤其是AIS數(shù)據(jù)，不僅數(shù)據(jù)量大，而且數(shù)據(jù)之間存在相互關(guān)聯(lián)。[1]結(jié)合數(shù)據(jù)挖掘技術(shù)，在海上交通工程的基礎(chǔ)上，將船舶位置空間、時間、船型和船速等屬性有機融合進行聯(lián)合聚類，以識別重點水域,發(fā)現(xiàn)典型特征，為實現(xiàn)智能水域監(jiān)管和海上安全保障提供理論基礎(chǔ)。

在對交通特征的研究中，時空數(shù)據(jù)的耦合特性逐漸受到關(guān)注。文獻[2]利用信息采集軟件和數(shù)據(jù)挖掘技術(shù)對大量的AIS數(shù)據(jù)進行分析，識別特定水域的船舶會遇態(tài)勢，并分析該水域船舶會遇空間分布和時間分布。文獻[3]引入k-means算法對數(shù)據(jù)設(shè)定時間維度上興趣時間段的約束條件，提取興趣位置點，實現(xiàn)對車與地理位置關(guān)系的挖掘。文獻[4]根據(jù)船舶交通特點將船舶領(lǐng)域概念引入DBSCAN算法中，對比船舶簇的速度和交通流的速度，識別出擁擠區(qū)域。文獻[5]結(jié)合AIS數(shù)據(jù)的具體特征提出時間切片化方法，在DBSCAN算法的基礎(chǔ)上，綜合考慮時間和空間要素，提出船載AIS數(shù)據(jù)時空聚類算法，并對實際數(shù)據(jù)進行分析。時空聚類現(xiàn)已成為海量時空數(shù)據(jù)分析的一個重要手段和前沿研究方向。[5-6]

船載AIS數(shù)據(jù)屬于典型的時空數(shù)據(jù)，其所包含的船舶空間、時間和其他維度屬性數(shù)據(jù)中蘊含著大量的潛在特征，單一維度的聚類很難發(fā)現(xiàn)數(shù)據(jù)中潛在的有用信息，多屬性聯(lián)合聚類分析已成為船載AIS數(shù)據(jù)挖掘的趨勢?，F(xiàn)有研究傾向于在已有聚類算法上，尤其是對DBSCAN算法進行改進，從而對海上交通特征的時空關(guān)系進行挖掘與分析，但當數(shù)據(jù)量增大時，聚類收斂時間變長，實現(xiàn)算法需消耗的內(nèi)存也變大，且參數(shù)選取相對復雜。由此，本文提出利用非負矩陣分解方法對船載AIS數(shù)據(jù)進行時空聯(lián)合聚類，挖掘海上交通特征和船舶行為模式，并聯(lián)系實際進行分析。

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

從海量數(shù)據(jù)中提取關(guān)注的信息，首先需對數(shù)據(jù)進行過濾和清洗。[7]在實際應用中，從數(shù)據(jù)庫中獲取的原始船載AIS數(shù)據(jù)集不能直接利用，需進行數(shù)據(jù)清洗與預處理。預處理的目的是將原始粗糙的數(shù)據(jù)轉(zhuǎn)化為適合分析處理的形式，其過程有特征選擇、維規(guī)約等。[8]

從AIS數(shù)據(jù)庫中提取指定區(qū)域和指定時間段內(nèi)的船舶靜態(tài)信息和動態(tài)信息。本文選取MMSI、經(jīng)度、緯度、UNIX時間戳、船型和航速等信息進行研究。

1) 船舶軌跡數(shù)據(jù)集T是指船舶水上航行所形成軌跡的數(shù)據(jù)信息集合，有

T={(time1,lon1,lat1),(time2,lon2,lat2),…,

(timeN,lonN,latN)}

(1)

式(1)中：timeN為第N個數(shù)據(jù)的UNIX時間戳；lonN為經(jīng)度；latN為緯度。

2) 時間域是指根據(jù)固定時間段將研究的時間切片化，并做編號標記，時間域為

t={t1,t2,…,tM}

(2)

式(2)中：tM為第M個時間域的編號。

3) 區(qū)域是指根據(jù)水域的形狀和航道，將研究區(qū)域劃分成的若干個不規(guī)則的小區(qū)域，有

R={r1,r2,…,rM}

(3)

式(3)中：rM為第M個區(qū)域的編號。

4) 標記數(shù)據(jù)集是指將原始數(shù)據(jù)集各屬性數(shù)據(jù)的值分別映射到區(qū)域、時間域上的標記組成的數(shù)據(jù)集合。

本文中提出一個船載AIS數(shù)據(jù)預處理框架，具體步驟為

1) 從AIS數(shù)據(jù)庫中提取指定時間與區(qū)域間的船舶動態(tài)信息、靜態(tài)信息和航次信息，形成初始軌跡數(shù)據(jù)集。

2) 篩除初始軌跡數(shù)據(jù)集中存在的錯誤、缺失的數(shù)據(jù)，將同一MMSI的AIS數(shù)據(jù)按時間戳排序并連接起來；以時間、航向和停靠點為依據(jù)，劃分出在同一MMSI軌跡數(shù)據(jù)中的不同航次和進出港情況。

3) 研究的水域劃分為若干個區(qū)域，將經(jīng)度和緯度數(shù)據(jù)映射到區(qū)域中并標記；研究的時間段劃分為若干個時間域，將UNIX時間戳轉(zhuǎn)換為標準時間，樣本數(shù)據(jù)映射到時間域后標記。

4) 構(gòu)建標記數(shù)據(jù)集。

2 稀疏約束下的非負矩陣分解

2.1 非負矩陣分解

矩陣分解是實現(xiàn)大規(guī)模數(shù)據(jù)處理與分析的一種有效方式[9]，通過矩陣分解解決實際問題的方法有很多，如主成分分析(Principal Component Analysis,PCA)、獨立成分分析(Independent Component Analysis,ICA)和奇異值分解(Singular Value Decomposition,SVD)等。[10]這些方法可通過迭代分解的方式，近似地將數(shù)據(jù)原始矩陣分解為2個較低維數(shù)矩陣的乘積，獲得數(shù)據(jù)的特征，各有優(yōu)勢，但存在一個共同、明顯的缺點，就是在分解過程中沒有對分解的結(jié)果進行非負約束限制，這些算法分解出的結(jié)果通常存在負數(shù)，從數(shù)學計算的角度看該結(jié)果是正確的，但在實際應用中，分解得到的負值結(jié)果通常是沒有意義的且無法解釋。LEE等[11]提出非負矩陣分解(Nonnegative MaxtrixFactor,NMF)，其是一種新的矩陣分解的思想。該算法通過對基與系數(shù)施加非負約束，保證分解結(jié)果均為正值；基矩陣利用系數(shù)矩陣作為權(quán)重，重構(gòu)原始矩陣。非負矩陣分解見圖1，其基本原理如下：

對于一個m×n的非負矩陣Vm×n，存在一個非負矩陣Wm×k和一個非負矩陣Hk×n，有

Vm×n≈Wm×kHk×n

(4)

即

Wmk≥0;Hkn≥0

(5)

圖1 非負矩陣分解示意

V中的列向量可解釋為對基矩陣W中所有列向量的加權(quán)和，而權(quán)重系數(shù)為系數(shù)矩陣H中對應列向量中的元素。非負矩陣分解對數(shù)據(jù)的學習有部分組成整體的特性，因為分解所得元素只有“純加性”，原始數(shù)據(jù)的整體由局部特征表示。由于其求解方法收斂速度快、左右非負矩陣存儲空間小，因此,其能實現(xiàn)高維的數(shù)據(jù)矩陣降維處理，適合處理大規(guī)模數(shù)據(jù)。目前，NMF算法在聚類/數(shù)據(jù)挖掘、特征學習、圖像分析、語音處理和生物醫(yī)學工程等方面有很好的應用，在海事領(lǐng)域已有SAR圖像船舶檢測[12-13]、潮汐數(shù)據(jù)的分析與預測[14]等方面的應用研究。

2.2 非負矩陣分解迭代算法

非負矩陣分解是求解出2個非負矩陣W∈Rm×r和H∈Rr×n,使重構(gòu)矩陣W×H與原始矩陣V之間的誤差最小，實際上是一個最優(yōu)化問題。問題的目標函數(shù)有很多，其中應用最廣泛的是LEE等[15]提出的基于歐氏距離和Kullback-Leibler離散度的2種度量模型，衡量重構(gòu)矩陣的近似化程度。

基于歐氏距離的目標函數(shù)定義為

(6)

基于Kullback-Leibler離散度的目標函數(shù)定義為

(7)

應用基于投影梯度法的交替最小二乘法[16]可求解非負矩陣分解問題。交替最小二乘法的算法模型如下：

1) 初始化W1,ia≥0,H1,bj≥0,i,a,b,j

2) Fork=1,2,…

(8)

(9)

該算法是在迭代中分別求解W和H,直至滿足終止條件。投影梯度法的迭代規(guī)則[17]為

崩塌在SPOT5的影像解譯特征：崩塌堆積體的平面形態(tài)多為弧形、扇形等。崩塌體后緣在遙感影像上陽坡為淺色調(diào)區(qū)塊、陰坡呈濃重的陰影區(qū)帶。崩塌堆積體上無植被覆蓋，見圖6(a)。

Hk+1=PΩ[Hk-ηHKPHK]

(10)

Wk+1=PΩ[Wk-ηWKPWK]

(11)

式(10)和式(11)中:PΩ[ξ]為ξ到凸集Ω={ξ∈R;ξ≥0}的一個映射;PHk和PWk分別為W和H的下降方向;ηHk和ηWk為相應的下降步長。該算法在迭代中同時求解W和H，直至滿足終止條件，可達到收斂速度更快的效果。

2.3 稀疏約束的非負矩陣分解

“稀疏編碼”旨在用少量的元素表示大量的有用數(shù)據(jù)，稀疏矩陣存儲信息更高效，占用的資源少，且利于解釋非負矩陣分解的結(jié)果。同類數(shù)據(jù)按列組成矩陣,矩陣的行向量對應某一特征,稀疏約束使其選擇盡量少的非零行向量達到特征選擇的目的。[18]非負矩陣分解得到的結(jié)果本身存在一定的稀疏性，但不夠充分，因此，在NMF中加上稀疏約束，可提高分解結(jié)果的質(zhì)量。本文通過在基于投影梯度法的交替最小二乘算法中，對基矩陣添加L0范數(shù)約束，提高NMF分解結(jié)果的稀疏度。問題的目標函數(shù)改寫為

(12)

(13)

式(13)中：參數(shù)L是Wi中允許的最大非零項數(shù)。

基于以上討論，提出在L0范數(shù)約束下，利用基于投影梯度法的交替最小二乘法的非負矩陣分解聚類算法為

Input：

標記數(shù)據(jù)集X∈Rm×n

結(jié)構(gòu)參數(shù)options：

分解得到矩陣的秩K

公差tolerance

外部循環(huán)次數(shù)numIter

最大迭代數(shù)maxIter

最大非零項數(shù)L

output：

分解后的基矩陣W∈Rm×K，系數(shù)矩陣H∈RK×n

Steps:

1) 隨機初始化矩陣H

2) 執(zhí)行外部循環(huán)：

利用非負最小二乘求解基于H對X的非負分解結(jié)果WT

fori，…，k

將Wi中最小的D-L個數(shù)值更新為零

根據(jù)得到的W，對H進行重新編碼

矩陣梯度=WT(WH-X)

(14)

根據(jù)式(14)進行迭代，直至滿足條件

3 實例分析

3.1 數(shù)據(jù)預處理

本文選取天津港水域真實AIS數(shù)據(jù)作為研究對象。

1) 研究海域為：經(jīng)度117°35′35″E ～ 118°43′E，緯度38°3′30″N ～ 39°34′30″N。參考交通流走勢，將該區(qū)域劃分為38個不規(guī)則區(qū)域，劃分詳情見圖2。

圖2 AIS樣本和區(qū)域劃分結(jié)果

2) 研究時間為：2018-01-01T00:00—2018-01-06T24:00?？紤]到研究海域的大小和船舶航速，在該時間范圍內(nèi)，按24 h制，每4 h劃分為一個時間域，共劃分36個時間域，編號詳細情況見表1。

表1 時間域劃分與編號結(jié)果 h

在經(jīng)過數(shù)據(jù)清洗之后，獲得998艘船舶的共 185 818條樣本數(shù)據(jù)，處理前的AIS數(shù)據(jù)列表見表2，處理后的標記數(shù)據(jù)列表見表3。

表2 處理前的AIS數(shù)據(jù)列表 (°)

表3 處理后的標記數(shù)據(jù)列表

3.2 時空數(shù)據(jù)聚類結(jié)果與分析

將AIS數(shù)據(jù)分為進港、出港和全部航次樣本數(shù)據(jù)集，分別進行分析。以時間-空間屬性矩陣為例,試驗采用SQLite 3.5作為后臺數(shù)據(jù)庫，利用PyCharm 2019進行數(shù)據(jù)清洗和預處理工作，采用MATLAB R2015 a 完成聚類算法的開發(fā)。船舶時空數(shù)據(jù)集未進行數(shù)據(jù)頻率上的約減，即同一艘船舶在時間段tM內(nèi)，在區(qū)域rM內(nèi)產(chǎn)生的所有符合研究條件的時空數(shù)據(jù)均納入研究數(shù)據(jù)集中，計入?yún)^(qū)域rM在時間段tM內(nèi)的船舶流量密度。

當L=4時，L影響分解結(jié)果的稀疏度和誤差，通過反復試驗，結(jié)合分解結(jié)果的可解釋性，選取L=4，基矩陣的稀疏度為0.690 7。當L=4時，對K的不同取值做30次獨立重復試驗，得到各K值下的平均誤差見表4。

由表4可知:當K值增大時，矩陣分解結(jié)果的誤差減小。對于此次試驗的數(shù)據(jù)集，當K≥6時，數(shù)據(jù)矩陣分解的誤差變化小于0.03。K值決定聯(lián)合聚類得到的共簇數(shù)量，當K值不斷增大時，可能會造成結(jié)果冗余。綜上所述，選取K=6、L=4(即原始數(shù)據(jù)集聚類得到的時空共簇數(shù)目為6，每個類中包含的屬性項維數(shù)為4)，分別對全部航次時空數(shù)據(jù)和出港船舶時空數(shù)據(jù)進行聚類分析。

表4 不同K值下的試驗平均誤差值

3.2.1以位置屬性為分析對象的聚類結(jié)果分析

為盡可能多地得到研究海域中的空間特征，使用全部航次的時空數(shù)據(jù)矩陣進行聚類分析，各類中的區(qū)域標記示意見圖3。

a)類簇1

由圖3可知：聚類結(jié)果結(jié)合實際的海域功能區(qū)可得到以下信息：

(1) 類簇1包括區(qū)域15、區(qū)域16、區(qū)域17和區(qū)域19，這是一條由區(qū)域19經(jīng)過大沽口北錨地進入天津港港口的一段航路；

(2) 類簇2包括區(qū)域10、區(qū)域18、區(qū)域20和區(qū)域30，此區(qū)域為天津港主航道出港航路；

(3) 類簇3包括區(qū)域9、區(qū)域19、區(qū)域23和區(qū)域25，此區(qū)域為天津港主航道進港航路；

(4) 類簇4包括區(qū)域6、區(qū)域8、區(qū)域22和區(qū)域24，這是一條進出曹妃甸港區(qū)的航路；

(5) 類簇5包括區(qū)域12、區(qū)域14、區(qū)域27和區(qū)域33，這是一條途徑區(qū)域12、區(qū)域33和大沽口航道進出大沽口港區(qū)的航路；

(6) 類簇6包括區(qū)域3、區(qū)域13、區(qū)域15和區(qū)域16，此區(qū)域為進出天津港東疆、南疆和北疆港區(qū)的航路。

圖3中，類簇1和類簇3中均包含區(qū)域19，類簇1和類簇6中均包含區(qū)域15和區(qū)域16，出現(xiàn)這種結(jié)果的原因有2個：

(1) 與區(qū)域的劃分有關(guān)，即1個區(qū)域中涵蓋2種船舶運動類型的區(qū)域，如區(qū)域15和區(qū)域16，進出港口的航道均在這2個區(qū)域內(nèi)；

(2) 與聚類數(shù)K和每一類中屬性項數(shù)L的取值有關(guān)。以類簇1和類簇3為例，2個類代表的區(qū)域見圖4。由圖4可知:類簇1和類簇3組成一段完整的航路，但在聚類表示中，由于選取的每一類中的屬性項數(shù)為4，即每一個類由最相關(guān)的4個區(qū)域表達，而該航路經(jīng)過的區(qū)域數(shù)量大于4，實際是通過類簇1和類簇3這2個類進行表示，因此，在對聚類結(jié)果上出現(xiàn)區(qū)域的重疊。

圖4 類簇1和類簇3涵蓋區(qū)域示意

從聚類的層面看，NMF屬于“軟聚類”算法，1個元素可屬于多種類型；從實際意義的角度分析，同一個區(qū)域中可能有一定比例的、不同運動模式的船舶航行。因此，在利用時空數(shù)據(jù)聚類對船舶與位置屬性的整體關(guān)系進行分析時，1個區(qū)域?qū)儆诙鄠€類是符合水上交通的實際情況的；從數(shù)值關(guān)系的層面看，在類簇1和類簇3中，區(qū)域19的船舶占比之和與船舶進入該區(qū)域的前段航路(區(qū)域23)占比接近，因此，該航路可看成被“拆分”在2個類簇中表達，此現(xiàn)象也符合NMF算法“由局部構(gòu)成整體”的特性。

3.2.2以時間屬性為分析對象的聚類分析

以研究船舶出港流量變化的時間規(guī)律為例，使用出港船舶的時空數(shù)據(jù)矩陣進行聚類分析，時空聚類結(jié)果的區(qū)域標記示意見圖5。

a)類簇1

以類簇3為例，對該類中的區(qū)域在時間模式下的波動情況進行分析，詳細情況見圖6。圖6中：每條線代表一個區(qū)域；x軸為1～6日劃分的36個時間域；y軸為時間段內(nèi)區(qū)域中的船舶流量數(shù)值；類簇3中包含的區(qū)域為區(qū)域8、區(qū)域10、區(qū)域22和區(qū)域30，包括天津港出港主航道的一段航路和進出曹妃甸港區(qū)的一段航路。

圖6 類簇3聚類關(guān)系中船舶流量隨時間波動的情況

由圖6可知:在類簇3包含的區(qū)域關(guān)系中，區(qū)域8、區(qū)域10、區(qū)域22和區(qū)域30的船舶出港流量變化趨勢接近，說明有一部分船舶航行途經(jīng)該類的區(qū)域8、區(qū)域10、區(qū)域22和區(qū)域30。結(jié)合實際分析，其反映的隱含關(guān)系是曹妃甸港區(qū)和附近水域船舶駛離曹妃甸時的航行方法，即由開敞性泊位區(qū)(區(qū)域8)或東側(cè)錨地(區(qū)域22)內(nèi)離泊，通過第三通航分道，駛向警戒區(qū)(區(qū)域30內(nèi))，由第一通航分道(區(qū)域10)內(nèi)駛離曹妃甸港區(qū)。由圖6可知：在1～6日這6 d中，船舶以該航行方式出港的峰值也呈現(xiàn)一定的規(guī)律，即在每天的8～12 h時間段內(nèi)均處于當日船舶出港最高峰，夜間航行出港的船舶數(shù)量相對較少。

4 有待進一步研究的問題

4.1 聚類特征數(shù)K的選取

目前，對于確定聚類特征數(shù)K以獲得最佳的聚類效果尚沒有很好的解決方法，因為,K的取值與數(shù)據(jù)集的類型有關(guān)，在對非負矩陣進行分解時，選取K值的約束條件為

K

(15)

式(15)中：當K較大時，降維效果不顯著，特征不突出；當K?min{m,n}時，特征明顯，但易導致忽略一些特征信息。目前,通用的方法是選取一系列K值分別進行試驗，將獲得的最佳辨識結(jié)果的聚類簇數(shù)作為K的取值。

4.2 區(qū)域劃分

對于區(qū)域劃分問題，特征聚類結(jié)果的可讀性受“規(guī)則網(wǎng)格區(qū)域”或“不規(guī)則區(qū)域”的劃分方式的影響不大，更多的是與區(qū)域劃分精度相關(guān)，例如在區(qū)域?qū)傩缘木垲惙治鲋?，區(qū)域劃分精度低，會出現(xiàn)1個區(qū)域分別隸屬于多個類簇的情況，例如圖3中區(qū)域15和區(qū)域16均隸屬于類簇1和類簇6，這種結(jié)果在規(guī)則和不規(guī)則的區(qū)域劃分中均會出現(xiàn)，但在低精度的區(qū)域劃分中出現(xiàn)的頻率更高；若劃分精度過高，矩陣特征維數(shù)增加，K值和每個類中特征維數(shù)L的選取難度增大。

4.3 最優(yōu)化問題

非負矩陣分解是一個約束優(yōu)化問題，矩陣分解的迭代收斂速度與目標函數(shù)的選取和稀疏約束的方式有關(guān)，因此,目標函數(shù)與稀疏約束函數(shù)等的優(yōu)化問題是進一步研究的方向。

5 結(jié)束語

海上交通時空數(shù)據(jù)存在耦合性，利用這種特性研究海上交通數(shù)據(jù)的時空分布特征和規(guī)律可發(fā)現(xiàn)其隱含的信息，對航道規(guī)劃和水上安全等工作的正常開展具有重要的現(xiàn)實意義。本文以計算機數(shù)據(jù)挖掘為研究手段，提出采用非負矩陣分解算法對大量的AIS數(shù)據(jù)在時間和空間2個屬性上同時進行約束聚類，探索海上交通數(shù)據(jù)時間和空間2個維度之間的關(guān)系，利用真實的AIS大數(shù)據(jù)集進行試驗，并對已有算法進行L0范數(shù)約束的改進，以達到更好的突出特征的效果，最終得到研究水域的主要航路，不同區(qū)域之間相近的船舶流量波動關(guān)系，以及針對途徑區(qū)域8、區(qū)域10、區(qū)域22和區(qū)域30駛離曹妃甸港區(qū)的船舶通常在8～12 h達到離港峰值等信息，聚類結(jié)果與實際相符，證明方法是可行的。該方法可用于挖掘研究水域的船舶行為模式，分析船舶的運動規(guī)律等，為水上交通安全監(jiān)管和海上安全保障相關(guān)研究提供一種新的思路。