蘇俊杰, 蘭培真
(1. 集美大學(xué)海上交通安全研究所, 福建 廈門 361021;2. 交通安全應(yīng)急信息技術(shù)國家工程實驗室, 福建 廈門 361021)
航運經(jīng)濟的發(fā)展促使船舶數(shù)量日益增加。船舶自動識別系統(tǒng)(automatic identification system,AIS)作為船與船、船與岸之間船舶相關(guān)信息的交換設(shè)備在船舶監(jiān)管中廣泛應(yīng)用,產(chǎn)生了海量的船舶動靜態(tài)信息數(shù)據(jù),而提取特征相似的船舶運動軌跡簇是分析船舶交通規(guī)律和船舶行為模式的基礎(chǔ)[1],因此,急需開展船舶軌跡的聚類研究。早期軌跡聚類研究[2-3]以軌跡點為樣本,忽略了相鄰軌跡點的時空關(guān)聯(lián),因此近期研究樣本大多為由連續(xù)點構(gòu)成的軌跡段。陳錦陽等[4]通過設(shè)定Hausdorff距離閾值劃分軌跡類別,類似指標還有歐氏距離、最長公共子序列距離、動態(tài)時間規(guī)整(dynamic time warping,DTW)距離等,然而這些指標只考慮船位偏差,描述船舶相似運動的能力有限。劉磊等[5]考慮船位、航速和航向定義綜合距離,并以K最鄰近分類器劃分軌跡;王維剛等[6]以8個軌跡特征度量軌跡相似性,設(shè)計了加權(quán)樸素貝葉斯分類器,證明了考慮多維特征度量相似性能更好地捕捉不同軌跡的差異[7]:這些監(jiān)督分類方法在航路固定時表現(xiàn)較好,但要求提供準確的先驗知識且難以發(fā)現(xiàn)隱蔽軌跡簇,無法滿足多數(shù)聚類情況的需要。因此,以具有噪聲的基于密度的空間聚類(density-based spatial clustering of applications with noise, DBSCAN)算法[8-11]為典型的無監(jiān)督聚類方法的應(yīng)用較多。然而,全局性的經(jīng)驗參數(shù)難以在密度不均的船舶交通環(huán)境下對各局部進行合理的聚類,為此,趙梁濱等[12]通過層次聚類獲取局部密度環(huán)境的適應(yīng)參數(shù)。此外,LI等[13]在K均值聚類算法的基礎(chǔ)上提出自適應(yīng)K值的多步軌跡聚類算法,但該算法直接處理整條軌跡,時空復(fù)雜度較高且丟失了軌跡的局部相似性;馬文耀等[14]采用譜聚類(spectral clustering,SC)算法區(qū)分船舶在不同行為模式下的軌跡;牟軍敏等[15]提出船舶軌跡的自適應(yīng)SC算法,但船位偏差仍是相似性度量的唯一依據(jù);冷泳林等[16]以線型特點衡量軌跡的相似性,并以仿射傳播(affinity propagation,AP)算法進行聚類,但聚類后的船舶運動軌跡相似性較低。針對上述船舶軌跡聚類方法存在的不足,本文提出一種考慮多維特征的船舶軌跡分層聚類(multi-features hierarchical clustering, MHC)算法,主要由分層軌跡相似性度量和CFA-DBSCAN算法2部分構(gòu)成,這里CFA(core firefly algorithm)指核心螢火蟲算法。
船舶運動慣性大,船舶行為模式具有區(qū)域性[10],相似軌跡應(yīng)為相似水域內(nèi)軌跡線型特點相似且具有相似航行動態(tài)的軌跡,因此以地理網(wǎng)格和時段作為軌跡段的劃分依據(jù)。任意軌跡可表示為T=(p1,p2,…,pN),其中:N為軌跡點數(shù);pi為i時刻的船舶運動狀態(tài),是由經(jīng)度λi、緯度φi、航速vi、航向αi和時隙Δti構(gòu)成的5維特征向量,即pi=(λi,φi,vi,αi,Δti)。軌跡T的線型特點及軌跡所處水域環(huán)境可由軌跡線型圖F1和水域信息圖F2描述,其中:F1為各軌跡點按順序連接而成的曲線圖;F2為軌跡所處地理網(wǎng)格范圍內(nèi)的地圖圖像。以圖1為例說明水域信息圖的構(gòu)建過程。為保證聚類結(jié)果的合理性,應(yīng)盡量選擇噪聲點少的,即同類地理環(huán)境區(qū)分度低、異類地理環(huán)境區(qū)分度高的地圖圖像數(shù)據(jù)。因為遙感影像包含的噪聲點較多,容易降低不同水域環(huán)境的辨識敏感度,所以本文研究區(qū)域的水域信息圖采用如圖1a所示的公開地圖圖像數(shù)據(jù)。地理網(wǎng)格的設(shè)置應(yīng)使得軌跡段的點數(shù)和水域信息圖采樣范圍合理,采用圖1b所示5×5地理網(wǎng)格將原始軌跡劃分到18個水域,即首先使由地理網(wǎng)格劃分的各軌跡段包含軌跡點數(shù)在[500,1 000]內(nèi),其次使地理網(wǎng)格內(nèi)的地圖圖像能達到粗略區(qū)分港口、沿岸和開闊水域的目的,則軌跡所處水域信息可用該水域范圍內(nèi)的地圖圖像表示。

a)研究區(qū)域的船舶軌跡

b)地理網(wǎng)格劃分的水域
將船舶的多維特征劃分為靜態(tài)圖像上層和動態(tài)數(shù)值下層,對不同層分別建立不同的軌跡相似性度量指標,經(jīng)兩次遞進聚類逐步放大軌跡簇間的差異,最終實現(xiàn)隱蔽軌跡簇的辨識。MHC基本框架見圖2,其中,SIFT(scale-invariant feature transform)指尺度不變特征變換算法。

圖2 MHC基本框架
上層聚類特征為軌跡線型圖和水域信息圖,軌跡線型圖和水域信息圖均為64×64的灰度圖像,各像素點取值區(qū)間為[0,255]。SIFT算法[17]魯棒性強,在光照、旋轉(zhuǎn)、視點和仿射變換條件下均能保持優(yōu)良的穩(wěn)定性,故采用SIFT算法建立上層相似性的度量指標。如圖3所示,SIFT算法本質(zhì)上是將圖像相似性轉(zhuǎn)換為2張圖像所有匹配特征點的相似性,因此,需要檢測圖像包含的所有特征點,并為每個特征點生成描述符。

圖3 SIFT算法基本框架
特征點檢測階段需要確定特征點所在的尺度和位置。以圖4為例說明圖像特征點的檢測過程。給定圖像I(x,y)和尺度高斯變換函數(shù)G(x,y,σ)(σ為采樣尺度),可構(gòu)建圖像的尺度空間L(x,y,σ)。尺度空間是由不同分辨率圖像構(gòu)成的高斯金字塔,用于描述圖像的尺度不變性。
(1)
L(x,y,σ)=G(x,y,σ)*I(x,y)
(2)
式中:σ取1.6;*為卷積運算。
高斯差分尺度空間D(x,y,σ)為相鄰尺度的高斯差分與原圖像的卷積,特征點由高斯差分尺度空間的局部極值點構(gòu)成。為尋找局部極值點,一個像素點要與其所有的相鄰點比較,看該像素點是否比其圖像域和尺度域的相鄰點大或者小。如圖4c所示,檢測點需要與26個點(與檢測點同尺度圖像域的8個相鄰點+與檢測點上下相鄰尺度圖像域的9×2個點)進行比較,以確保在尺度空間和圖像空間都檢測到極值點。
D(x,y,σ)=
(G(x,y,kσ)-G(x,y,σ))*I(x,y)
(3)

a)高斯金字塔

b)高斯差分金字塔

c)特征點
特征點描述階段需要為每個特征點生成對應(yīng)的特征描述符,即以特征點鄰域像素的梯度方向分布賦予特征點方向信息,L尺度下(x,y)處的梯度值m(x,y)和方向θ(x,y)分別為
m(x,y)=((L(x+1,y)-L(x-1,y))2+
(L(x,y+1)-L(x,y-1))2)1/2
(4)
(5)
以圖5為例說明特征描述符在實際應(yīng)用中的構(gòu)建。首先將特征點的鄰域像素劃分為2×2的區(qū)塊,每個區(qū)塊由4×4子域構(gòu)成,計算每個子域的梯度并統(tǒng)計各區(qū)塊內(nèi)8個方向(每個方向的范圍為0~45°)的梯度方向直方圖,其中,縱軸為每個梯度方向的梯度累計值,則每個區(qū)塊可描述為一個8維的向量,該特征點的特征描述符為一個2×2×8維的向量。

a)圖像梯度

b)統(tǒng)計直方圖

c)關(guān)鍵點描述
特征點匹配旨在尋找兩張圖像匹配的特征點對,匹配度量為特征描述符的余弦距離。匹配特征點對的查找方法:檢索某圖像特征點在另一圖像中的最鄰近和次鄰近特征點,若最近距離與次近距離的比值小于給定閾值(通常取0.6~0.8[18]),則該特征點與最鄰近特征點匹配。圖像距離為所有匹配特征點對余弦距離的平均值。上層船舶軌跡相似性可用軌跡線型圖和水域信息圖圖像距離的均值Dimg度量。
(6)
式中:d1為軌跡線型圖的圖像距離;d2為水域信息圖的圖像距離。
下層聚類特征為船位、航速、航向和時隙,軌跡相似性即軌跡運動特征的相似性。船舶軌跡是不等長的時序數(shù)據(jù),因此,采用DTW[10]算法建立下層軌跡相似性的度量指標。以圖6為例說明下層軌跡相似性的計算過程:給定任意兩條船舶軌跡TA和TB,構(gòu)建初始距離矩陣,該矩陣的元素d(i,j)為對應(yīng)兩軌跡點特征向量的余弦距離。余弦距離按式(7)計算,計算前需將各特征分量進行歸一化以消除不同量綱的影響。按式(8)計算初始距離矩陣的累積距離矩陣,即累積距離矩陣中Mc(i,j)的值等于其左、上和左上位置數(shù)值中的最小值與初始距離矩陣對應(yīng)位置之和。下層軌跡的相似性可由遍歷軌跡TA和TB所有軌跡點的最小累積距離Mc衡量。
(7)
Mc(i,j)=d(i,j)+min(Mc(i-1,j-1),
Mc(i,j-1),Mc(i-1,j))
(8)
式中:xik和xjk分別為點pi和pj的第k個特征分量。通常d(i,j)越小,pi與pj越相似。假設(shè)pi=(1,1,0,0,1),pj=(1,1,1,0,1),則根據(jù)式(7)有

圖6 累積距離矩陣計算
DBSCAN算法是基于密度的無監(jiān)督聚類算法[19],辨識隱蔽簇的能力強且無須設(shè)定簇的數(shù)量,適用于船舶軌跡數(shù)據(jù)噪聲點多、船舶行為模式不確定性高的情況,相關(guān)定義如下:
定義1核心樣本:對某樣本li,若在給定距離半徑r的鄰域ε(li)中樣本的數(shù)量不低于給定閾值δ,則該樣本為核心樣本。
定義2密度直達:某核心樣本可密度直達其鄰域包含的所有樣本。
定義3密度可達:若li密度直達lk,lk密度直達lj,則li密度可達lj,即若某核心樣本li密度直達的樣本lk為核心樣本,則樣本li密度可達lk鄰域的所有樣本,密度可達具有傳遞性。
定義4密度相連:若核心樣本lk密度可達li和lj,則li與lj密度相連,即由密度可達的所有核心樣本的鄰域構(gòu)成的區(qū)域內(nèi)的樣本都密度相連。
DBSCAN算法旨在將密度相連的樣本聚成同簇,而在實際聚類過程中,樣本鄰域通常存在交集。為避免DBSCAN的重復(fù)查詢,僅聚類核心樣本,對同簇的核心樣本集對應(yīng)的鄰域樣本集取并集來進行簇的擴充。
改進的DBSCAN算法同樣對參數(shù)的設(shè)置比較敏感,因此,基于螢火蟲算法(firefly algorithm,F(xiàn)A)[20]來優(yōu)化參數(shù)。FA是新興的群智能優(yōu)化算法,概念簡單、優(yōu)化性能好,在參數(shù)調(diào)優(yōu)問題上應(yīng)用廣泛。FA基本步驟[21]如下:
步驟1初始化螢火蟲數(shù)量N,最大吸引度β0,光吸收系數(shù)γ,步長因子ρ,ρ∈[0,1],最大迭代次數(shù)tmax和搜索精度ξ。
步驟2隨機初始化螢火蟲位置,計算目標函數(shù)值S作為螢火蟲各自的最大熒光亮度I0。
步驟3計算群體中螢火蟲的相對亮度I和吸引度β,相對亮度決定螢火蟲的移動方向。
(9)
式中:rij表示螢火蟲i與j的距離;β0為最大吸引度,即rij=0時的吸引度。
步驟4更新螢火蟲被吸引后的位置。
xi(t+1)=xi(t)+β(xj-xi)+ρ(μ-0.5)
(10)
式中:t為迭代次數(shù);μ為[0,1]內(nèi)均勻分布的隨機因子。
步驟5更新螢火蟲的亮度。
步驟6若達到最大迭代次數(shù)或滿足搜索精度,則轉(zhuǎn)下一步;否則重復(fù)步驟3~6。
步驟7輸出全局極值點和個體最優(yōu)亮度值。
FA缺乏全局隨機性,易陷入局部最優(yōu)[20]且螢火蟲群位置的頻繁更新使得算法時空復(fù)雜度較高,因此,提出CFA。區(qū)別于FA,CFA只選取亮度排在前20%的螢火蟲作為核心螢火蟲進行位置更新,并隨機分配其余螢火蟲的位置來確保全局隨機性;在沒有局部更亮的個體時,核心螢火蟲會向種群的全局最優(yōu)點移動,最優(yōu)位置即核心螢火蟲的收斂位置。為避免隨著迭代次數(shù)增加,核心螢火蟲圍繞最優(yōu)位置振蕩,按式(11)分段設(shè)置步長因子。
(11)
式中:ρ0為初始步長因子,ρ0∈[0,1];lij為螢火蟲i與j在不同維度上的距離。
CFA-DBSCAN算法為改進DBSCAN算法與CFA參數(shù)尋優(yōu)算法的結(jié)合,即采用CFA優(yōu)化改進DBSCAN算法的參數(shù)r和δ以獲取當前最優(yōu)的軌跡聚類效果。為尋求參數(shù)r和δ的最優(yōu)值,定義參數(shù)優(yōu)化的目標函數(shù)為軌跡聚類的總輪廓系數(shù)S[22],即使得同簇內(nèi)樣本相似度最大,不同簇間的相似度最小[23]。其中,S取值越接近1說明算法聚類性能越好。
式中:s(i)為軌跡的輪廓系數(shù);a(i)為樣本i與同簇內(nèi)其他樣本的平均距離;b(i)為樣本i與不同樣本的平均距離;NS為樣本總數(shù)。
以廈門港及其附近水域的AIS數(shù)據(jù)來驗證MHC算法的有效性。AIS數(shù)據(jù)集包含船舶軌跡2 364條。為確保檢驗分析結(jié)果的準確性,參考ZHAO等[24]提供的方法對采集的AIS軌跡數(shù)據(jù)進行預(yù)處理,先用地理網(wǎng)格劃分船舶軌跡,再按小時時段劃分船舶軌跡,得到有效軌跡段共5 006條。該軌跡數(shù)據(jù)集的時間跨度為20190815T114425—20190820T114425,地理跨度為東經(jīng)117°54′37″~118°36′39″,北緯24°06′38″~24°31′37″,航速范圍為0~14 kn,航向范圍為0°~359.9°。
為檢驗MHC算法的性能,以異常軌跡數(shù)、簇數(shù)、簇內(nèi)DTW距離均值、簇間DTW距離均值和準確率為評價指標,對比分析MHC算法與CFA-DBSCAN、DBSCAN[25]、AP[16]、SC[14]、K均值聚類[13]等5種其他常用船舶軌跡聚類算法的評價指標值,見表1。5種其他常用算法均以船舶軌跡的DTW距離作為算法的樣本距離;CFA-DBSCAN算法的參數(shù)按螢火蟲種群規(guī)模N=100,最大迭代次數(shù)tmax=200,初始步長因子ρ0=0.8,光吸收系數(shù)γ=1設(shè)置;設(shè)定SC和K均值聚類算法的聚類簇數(shù)與MHC算法的聚類簇數(shù)相同。由表1可見,在無須事先設(shè)定簇數(shù)的聚類算法(MHC、DBSCAN和CFA-DBSCAN)中,MHC算法將軌跡聚為9簇,而DBSCAN算法和CFA-DBSCAN算法均只將軌跡聚為3簇,表明MHC算法能更好地辨識隱蔽軌跡簇。在將軌跡聚為3簇的算法中,CFA-DBSCAN算法的各指標均優(yōu)于AP算法和DBSCAN算法,這表明根據(jù)經(jīng)驗設(shè)置的DBSCAN算法參數(shù)不足以區(qū)分復(fù)雜的船舶軌跡,而通過CFA確定的最優(yōu)參數(shù)能在一定程度上提高DBSCAN算法的性能。在將軌跡聚為9簇的算法中,MHC算法在各評價指標上都優(yōu)于SC算法和K均值聚類算法,簇內(nèi)DTW距離均值為5.199,簇間DTW距離均值為18.032,準確率為91.50%,表明MHC算法能更好地辨識軌跡間的異同,從而實現(xiàn)船舶軌跡的準確聚類。

表1 各軌跡聚類算法的評價指標值
為分析MHC算法各模塊在軌跡聚類過程中的影響,設(shè)置3組并行實驗:第一組作為對照實驗,僅考慮船位、航速和航向特征,并用DBSCAN算法對軌跡進行聚類;第二組實驗將DBSCAN算法替換為CFA-DBSCAN算法;第三組實驗在第二組實驗的基礎(chǔ)上引入水域環(huán)境、軌跡線型和時隙特征對船舶軌跡進行分層聚類。如圖7所示:DBSCAN算法將船舶軌跡劃分為3簇,基本對應(yīng)著港區(qū)和外海的船舶,但局部軌跡分類明顯不符合實際,如外海船舶被歸到港區(qū)的船舶簇,且異常軌跡較多,聚類效果較差。這表明在較大水域環(huán)境范圍內(nèi)船舶軌跡通常密度不均,DBSCAN算法全局性經(jīng)驗參數(shù)的設(shè)置難以獲取局部密度環(huán)境下的合理聚類,且在區(qū)分船舶行為多變的軌跡方面存在局限性。CFA-DBSCAN算法同樣將船舶軌跡劃分為3簇,但異常軌跡數(shù)明顯減少,表明DBSCAN算法對參數(shù)的選取敏感,而通過CFA確定的參數(shù)能在一定程度上提高算法的質(zhì)量。MHC算法通過上層聚類后得到3簇軌跡,基本對應(yīng)著沿岸和外海曲度不同的船舶軌跡,異常軌跡相對較少,異常軌跡基本集中在水上交通環(huán)境復(fù)雜的港區(qū);通過下層聚類后,船舶軌跡被劃分為9簇,基本對應(yīng)船舶進出港區(qū)的軌跡匯聚和軌跡分流,與目前廈門的船舶交通狀況基本吻合,且異常軌跡明顯減少。這表明MHC算法通過上層和下層遞進聚類的方式能更好獲取局部密度環(huán)境的適應(yīng)參數(shù)。
表2為算法的運行指標。由表2可見:CFA-DBSCAN算法的鄰域查詢次數(shù)比DBSCAN的少,表明鄰域查詢部分的改進減少了聚類過程中的重復(fù)操作;CFA-DBSCAN算法聚類時間比DBSCAN算法的多,這是由于DBSCAN算法參數(shù)是事先設(shè)定的,不需要通過多次迭代更新,節(jié)省了聚類時間,而CFA-DBSCAN算法需要通過多次迭代使得軌跡聚類結(jié)果的輪廓系數(shù)收斂于全局最優(yōu);MHC算法需要的聚類時間比前兩種算法的多,這主要是因為在傳統(tǒng)船舶軌跡相似性考慮標準的基礎(chǔ)上引入了水域環(huán)境、軌跡線型和時隙特征,增加了軌跡相似性度量復(fù)雜度,而圖像相似度的計算和上下層遞進聚類增加了軌跡相似性度量的時間和鄰域查詢次數(shù),但也正是這些原因使得MHC算法更容易發(fā)現(xiàn)隱蔽軌跡簇,在聚類結(jié)果和準確度方面較優(yōu)。

a)DBSCAN

b)CFA-DBSCAN

c)MHC上層聚類

d)MHC下層聚類

表2 軌跡聚類算法的運行指標
為準確聚類復(fù)雜的船舶軌跡和辨識隱蔽軌跡簇,本文提出考慮多維特征的船舶軌跡分層聚類(MHC)算法。在傳統(tǒng)船舶軌跡聚類特征的基礎(chǔ)上引入水域環(huán)境、軌跡線型和時隙特征分層建立軌跡相似性度量指標并以CFA-DBSCAN算法遞進聚類各層軌跡,以廈門港及其附近水域AIS數(shù)據(jù)驗證了MHC算法的有效性。檢驗結(jié)果表明,MHC算法相較于常用的聚類算法能更準確地區(qū)分船舶軌跡并發(fā)現(xiàn)隱蔽軌跡簇,聚類后同簇內(nèi)軌跡的船舶運動特性更相似,而不同簇軌跡的船舶運動特性差異更明顯。聚類結(jié)果可為船舶交通流特性分析和船舶行為模式識別提供精確的樣本支持。MHC算法的聚類時間較長,因此,如何兼顧MHC算法聚類性能并優(yōu)化聚類時間將是后續(xù)研究的重點。