劉海楊,孟令航,林仲航,谷源濤
(清華大學電子工程系,北京 100084)
通過對海量歷史軌跡數據進行深度挖掘,發現未掌握的飛行器航路,進而掌握飛行器的運動規律,是統籌規劃全局空域交通管制的基礎,同時也是異常軌跡預警的關鍵。
學術界關于航路發現的研究大多是針對浮動車輛的GPS(Global Positioning System)軌跡數據,目前,未搜集到航路發現相關的研究。
基于浮動車輛GPS 軌跡數據的路網提取方法可以分為4類:
1)聚類方法。對軌跡點或線進行聚類,如Edelkamp等[1]、Schroedl等[2]利用K-means算法對軌跡進行聚類,采用樣條曲線擬合道路中心線來提取路網;Worrall 等[3]對軌跡點進行聚類,同時考慮方向角因素,鏈接各個節點提取路網;米春蕾等[4]、王冬等[5]對數據進行網格化后采用聚類方法,通過鏈接骨架節點或擬合中心線的方式,得到路網幾何數據。
2)Delaunay 三角網法。楊偉等[6]、唐爐亮等[7]對軌跡線構建約束三角網來提取路網。
3)基于柵格化和核密度估計的圖像處理方法。Shi等[8]、李俊杰等[9]對軌跡數據進行柵格化處理;Biagioni 等[10]采用核密度估計將軌跡數據轉換為圖像,然后采用圖像處理[11]相關手段提取路網。
4)增量合并方法。Zhang 等[12]、Bruntrup 等[13]、Cao 等[14]主要是根據已有路網和新的軌跡數據對路網進行更新。
飛行器飛行軌跡數據與浮動車輛的軌跡數據有著明顯差異:浮動車輛均按照現實存在的道路行駛,行駛線路較簡單;而飛行器的航路是抽象的,且飛行軌跡較復雜,前述基于浮動車輛GPS 軌跡數據的路網提取研究成果各有優點,但無法充分利用飛行器運動軌跡的特點和解決軌跡采樣率不同、高噪聲的問題,直接應用于飛行器航路發現時性能不佳。存在的主要問題如下:
1)雖然眾多研究成果對“直線”型航路段的提取效果較好,但對于“圓”狀航路結構和多分支的交叉路口結構的提取效果欠佳。“圓”狀航路結構容易提取為“點”結構,多分支結構中的部分夾角較小的分支容易被合并。
2)點、線聚類方法對采樣率和噪聲要求比較嚴格。
3)Delaunay 三角網、柵格化與核密度估計方法提取航路時,道路段呈鋸齒狀。
4)增量合并方法主要用于對航路的更新,而且對噪聲敏感。
因此,本文設計了一種基于軌跡點聚類的航路發現方法,旨在解決飛行器軌跡數據采樣率不同、高噪聲問題的同時,改善對特殊航路結構的提取效果。
本文設計了一種抗噪魯棒、易于實現的基于聚類的航路數據提取方法,具體步驟為:
1)對軌跡數據進行預處理,包括離群點剔除和卡爾曼濾波;
2)對預處理后的軌跡數據進行重采樣,使軌跡數據的采樣密度基本一致;
3)對軌跡數據點進行聚類,得到航路節點(聚類中心);
4)對航路節點進行修正,降低航路節點的冗余度;
5)根據軌跡數據點的連接關系,將航路節點進行連接,形成航路幾何數據。
本文涉及的軌跡數據由大量離散軌跡數據點組成,并存儲在數據庫中。軌跡數據點包括軌跡序號、經度、緯度、時間4 項有效屬性。由于各種因素的影響,在數據庫中存儲的軌跡數據存在噪聲和離群點,導致數據無法直接使用。針對軌跡數據中存在的噪聲和離群點,本文提出了如下解決方案。
1.1.1 離群點剔除
采用一種自適應閾值法對離群點進行剔除操作。由經度、緯度和時間計算軌跡點的平均速度,采用基于數理統計中的3σ 原則,對平均速度大于閾值的軌跡點進行剔除。
1)閾值初始化。根據式(1)計算出整條軌跡各個采樣點的平均速度。由于中值對噪聲不敏感,所以本文以平均速度的中值乘以一個增益系數η來初始化閾值。

式中:vi是pi點的速度,E0為初始化閾值,(xi,yi)為軌跡中第i個軌跡點pi的坐標,ti為軌跡點pi對應的時間。
2)由于當前軌跡點的運動狀態與相鄰軌跡點的運動狀態有一定的相關性,本文采用滑窗法對整條軌跡的采樣點進行判斷和剔除。該部分分為兩階段:
①窗未滿時,以初始化閾值為判斷依據。

式中:pi為軌跡的采樣點,當pi=正常點時,pi進入滑窗;否則,將pi從軌跡中刪除。
②窗滿時,根據3σ 原則,對窗內軌跡點的平均速度計算閾值。

式中:E為閾值,meanwin為窗內平均速度,varwin為窗內速度標準差。
當vi>E,pi=離群點,進行剔除;否則,窗體滑動,將窗體中最舊的點移除窗體,令當前點進入窗體,直至整條軌跡結束。
1.1.2 濾噪平滑
采用卡爾曼濾波對軌跡數據進行平滑處理,達到消除軌跡數據中部分誤差、平滑軌跡的效果。
1.2.1 孤立點剔除
Van Winden 等[15]通過統計分析得出結論:同一道路上的軌跡點到道路中心線的距離近似遵循高斯分布。軌跡預處理僅降低了軌跡數據的噪聲影響,但噪聲仍存在。由于部分軌跡數據中會存在孤立點,偏離航空主干道路的情況,為進一步降低對航路節點提取的影響,本文采用DBSCAN(Density-Based Spatial Clustering of Applications with Noise)算法[16]對軌跡數據點進行深度去噪處理。該算法不僅能夠對軌跡數據點進行聚類,且對噪聲不敏感,能處理非凸數據,較好地識別孤立點。
定義 孤立點。距離航路中心線較遠且稀疏的點,如圖1 所示。

圖1 孤立點示意圖Fig.1 Schematic diagram of isolated point
設軌跡{pi,1,pi,2,…,pi,n},其中pi,j為第i條軌跡的第j個軌跡點,軌跡集合I={T1,T2,…,Tm},軌跡點集合S={p1,1,p1,2,…,pm,n}。
1)對軌跡點集合S進行DBSCAN 聚類并設置合適的eps和Minpts,可得到聚類結果S1,S2,…,Sk,S-1,其中,S-1為孤立點集合,其余聚類結果視為有效點。
2)將孤立點集合S-1中的軌跡點從軌跡集合I中的軌跡T中剔除。
1.2.2 軌跡重采樣
在歷史軌跡數據庫中,軌跡數據的采樣時間間隔不統一,且同一軌跡的采樣點的距離間隔也不同,所以本文采用軌跡重采樣的方法,能消除采樣時間間隔不統一造成的軌跡點的空間密度不同帶來的影響,保證各條軌跡數據的軌跡點密度基本一致。該步驟的實現是為了1.2.3 節中的軌跡點聚類做好數據基礎。具體算法如下:


1.2.3 軌跡點聚類
對數據庫中的歷史數據進行處理后得到的每條軌跡大致為等距離采樣點,換言之,所有的軌跡數據點的密度較均勻,本文采用mini batchK-means 聚類對所有的軌跡采樣點進行聚類,根據數據規模和覆蓋面積,設定合適的聚類中心個數,經過聚類后,聚類中心能夠較好地反映出軌跡的航路結構,得到航路節點G={C1,C2,…,Cn}。
對航路節點進行提取后,由于上述聚類方法選取聚類中心個數的原因,會存在航路節點分裂現象,即兩個或多個節點之間的距離較近,本文采用基于密度的聚類方法,對航路節點進行合并,得到新的航路節點集合G′={C′1,C′2,…,C′n′},并對軌跡數據點pi所屬的節點簇類進行相應的修改,如圖2所示。

圖2 航路節點修正示意圖Fig.2 Schematic diagram of route node amendment
根據1.3 節可得到航路節點,將航路節點視為聚類中心,那么每條軌跡Ti={pi,1,pi,2,…,pi,n}的軌跡點pi,j(1≤j≤n)都含有各自的所屬類別,本文將軌跡Ti映射到航路節點集合中,每條軌跡就可由Ti={pi,1,pi,2,…,pi,n}轉換為由航路節點集合中的元素構成的序列為Ti={C′k1,C′k2,…,C′kn},其中:{pi,1,pi,2,…,pi,m}∈C′k1,{pi,m+1,pi,m+2,…,pi,m+l}∈C′ki,{pi,m+l+1,pi,m+l+2,…,pi,n}∈C′kn。
軌跡點pi映射到航路節點集合的原則:

式中:distance(p,C)為軌跡點p和航路節點C的歐氏距離。
將所有歷史軌跡映射為航路節點集合中的元素序列后,遍歷數據庫中的所有軌跡,可以得到整個航路節點的連接關系,進而得到航路幾何數據。
本文根據真實軌跡的運動特征和數據特點,生成仿真數據用于實驗。在生成仿真數據時,進行添加噪聲、離群點、軌跡斷裂等處理,使仿真數據盡可能逼近真實數據。實驗數據包括不同噪聲強度的軌跡數據,每種噪聲強度下包括4 800條軌跡。實驗數據集共有2 296 660 個軌跡數據點,每個軌跡點包括經度、緯度、時間等屬性。實驗數據集如圖3 所示。

圖3 實驗數據集Fig.3 Experimental data set
對實驗數據集進行預處理(如圖4 所示)后,通過對孤立點剔除(如圖5 所示),軌跡重采樣、軌跡點聚類和對航路節點的修正與連接,得到航路集合數據。針對實驗數據集,本文設定聚類中心個數為300,為了對航路節點的提取效果進行評價,本文以仿真生成實驗軌跡數據使用的航路結構作為參考,將本文的實驗結果與參考航路進行疊加對比分析。

圖4 軌跡預處理前后對比Fig.4 Comparison of trajectories before and after preprocessing

圖5 孤立點剔除效果Fig.5 Visualization of isolated point removal
2.2.1 定性評價
圖6(a)為本文方法的疊加分析圖。從整體結果來看,本文方法提取的航路節點基本覆蓋了參考航路,并且準確度較高;但對局部結果觀察發現,連接航路節點的過程中,部分航路會出現錯誤連接的情況,這主要是由于軌跡數據中存在的噪聲的影響,使得軌跡數據點向航路節點集合中映射出現偏差,導致錯誤連接。

圖6 兩種方法航路發現結果對比Fig.6 Comparison of results of two route discovery methods
2.2.2 定量評價
為了定量評價本文方法的實驗結果,將本文實驗結果與柵格化方法[8]的實驗結果分別與參考航路進行匹配,然后引用文獻[12]提出的緩沖區檢測方法,分別建立10 km、20 km、30 km 的緩沖區,對兩種方法提取的航路節點的覆蓋率α和航路長度的覆蓋率β進行性能評價。具體評價函數如下所示:


式中:pi為航路節點,N為航路節點總數,linei為提取的航路中第i段軌跡段在參考航路緩沖區的部分,LINEi為提取的航路中第i段軌跡段。
如圖7 和圖8 所示,本文方法提取的航路數據在節點覆蓋率和長度覆蓋率上有較大提高,特別是在10 km 的緩沖區內的實驗結果。對單一路段的提取,兩種方法的提取效果比較相近,但對于兩條航路夾角較小的交匯路口處,柵格化方法易提取成為一條航路;對“圓”狀航路結果容易誤提取為點結構。相較于柵格化方法,本文方法能夠在保持提取平滑航路準確率的基礎上,更好地提取“圓”狀航路結構和夾角較小的兩條航路的航路結構,有較好的抗噪性。

圖7 本文方法與柵格化方法的節點覆蓋率對比Fig.7 Comparison of node coverage between two methods

圖8 本文方法與柵格化方法的長度覆蓋率對比Fig.8 Comparison of length coverage between two methods
本文方法的有效性已在真實數據集上得到驗證,由于真實數據不便公開引用,所以本文對22 類民航軌跡數據進行航路提取來驗證本文方法的有效性。軌跡數據可視化如圖9(a)所示,航路提取結果如圖9(b)所示。根據民航數據的航路發現結果可知,本文方法能夠較好地發現航路,但存在一定的局限性,當某條航路或航路段上的軌跡數目稀少時,該航路或航路段上的軌跡易被檢測為孤立點進行剔除,因此漏掉部分航路。

圖9 民航軌跡數據實驗結果Fig.9 result of civil aviation data set
為了快速準確地從航空軌跡數據中提取航路的幾何數據,本文提出了一種基于軌跡點聚類的航路提取方法,利用航空軌跡模擬數據進行實驗分析,以模擬數據的理想航路為參考數據,與柵格化方法進行定性定量分析評價,結果表明本文方法提取的航路數據在節點覆蓋率和航路長度覆蓋率上都有明顯提高,而且本文方法易于實現,有一定的有效性和可靠性。
本文方法提取的航路節點,可將軌跡數據映射為節點序列,進而可將軌跡預測問題轉化為序列到序列的序列預測問題。為提高預測問題的正確率,需進一步研究和完善:1)軌跡數據映射到節點集合中的更有效的方法;2)各個航空路段的軌跡數據的密度相差較懸殊的情況,本文方法還需要進一步改進;3)由于噪聲的影響,連接航路節點時會有誤連的情況,下一步需對連接節點的方法進行進一步研究。