劉海楊，孟令航，林仲航，谷源濤

（清華大學(xué)電子工程系，北京 100084）

0 引言

通過(guò)對(duì)海量歷史軌跡數(shù)據(jù)進(jìn)行深度挖掘，發(fā)現(xiàn)未掌握的飛行器航路，進(jìn)而掌握飛行器的運(yùn)動(dòng)規(guī)律，是統(tǒng)籌規(guī)劃全局空域交通管制的基礎(chǔ)，同時(shí)也是異常軌跡預(yù)警的關(guān)鍵。

學(xué)術(shù)界關(guān)于航路發(fā)現(xiàn)的研究大多是針對(duì)浮動(dòng)車輛的GPS（Global Positioning System）軌跡數(shù)據(jù)，目前，未搜集到航路發(fā)現(xiàn)相關(guān)的研究。

基于浮動(dòng)車輛GPS 軌跡數(shù)據(jù)的路網(wǎng)提取方法可以分為4類：

1）聚類方法。對(duì)軌跡點(diǎn)或線進(jìn)行聚類，如Edelkamp等［1］、Schroedl等［2］利用K-means算法對(duì)軌跡進(jìn)行聚類，采用樣條曲線擬合道路中心線來(lái)提取路網(wǎng)；Worrall 等［3］對(duì)軌跡點(diǎn)進(jìn)行聚類，同時(shí)考慮方向角因素，鏈接各個(gè)節(jié)點(diǎn)提取路網(wǎng)；米春蕾等［4］、王冬等［5］對(duì)數(shù)據(jù)進(jìn)行網(wǎng)格化后采用聚類方法，通過(guò)鏈接骨架節(jié)點(diǎn)或擬合中心線的方式，得到路網(wǎng)幾何數(shù)據(jù)。

2）Delaunay 三角網(wǎng)法。楊偉等［6］、唐爐亮等［7］對(duì)軌跡線構(gòu)建約束三角網(wǎng)來(lái)提取路網(wǎng)。

3）基于柵格化和核密度估計(jì)的圖像處理方法。Shi等［8］、李俊杰等［9］對(duì)軌跡數(shù)據(jù)進(jìn)行柵格化處理；Biagioni 等［10］采用核密度估計(jì)將軌跡數(shù)據(jù)轉(zhuǎn)換為圖像，然后采用圖像處理［11］相關(guān)手段提取路網(wǎng)。

4）增量合并方法。Zhang 等［12］、Bruntrup 等［13］、Cao 等［14］主要是根據(jù)已有路網(wǎng)和新的軌跡數(shù)據(jù)對(duì)路網(wǎng)進(jìn)行更新。

飛行器飛行軌跡數(shù)據(jù)與浮動(dòng)車輛的軌跡數(shù)據(jù)有著明顯差異：浮動(dòng)車輛均按照現(xiàn)實(shí)存在的道路行駛，行駛線路較簡(jiǎn)單；而飛行器的航路是抽象的，且飛行軌跡較復(fù)雜，前述基于浮動(dòng)車輛GPS 軌跡數(shù)據(jù)的路網(wǎng)提取研究成果各有優(yōu)點(diǎn)，但無(wú)法充分利用飛行器運(yùn)動(dòng)軌跡的特點(diǎn)和解決軌跡采樣率不同、高噪聲的問(wèn)題，直接應(yīng)用于飛行器航路發(fā)現(xiàn)時(shí)性能不佳。存在的主要問(wèn)題如下：

1）雖然眾多研究成果對(duì)“直線”型航路段的提取效果較好，但對(duì)于“圓”狀航路結(jié)構(gòu)和多分支的交叉路口結(jié)構(gòu)的提取效果欠佳?！皥A”狀航路結(jié)構(gòu)容易提取為“點(diǎn)”結(jié)構(gòu)，多分支結(jié)構(gòu)中的部分夾角較小的分支容易被合并。

2）點(diǎn)、線聚類方法對(duì)采樣率和噪聲要求比較嚴(yán)格。

3）Delaunay 三角網(wǎng)、柵格化與核密度估計(jì)方法提取航路時(shí)，道路段呈鋸齒狀。

4）增量合并方法主要用于對(duì)航路的更新，而且對(duì)噪聲敏感。

因此，本文設(shè)計(jì)了一種基于軌跡點(diǎn)聚類的航路發(fā)現(xiàn)方法，旨在解決飛行器軌跡數(shù)據(jù)采樣率不同、高噪聲問(wèn)題的同時(shí)，改善對(duì)特殊航路結(jié)構(gòu)的提取效果。

1 技術(shù)路線與主要方法

本文設(shè)計(jì)了一種抗噪魯棒、易于實(shí)現(xiàn)的基于聚類的航路數(shù)據(jù)提取方法，具體步驟為：

1）對(duì)軌跡數(shù)據(jù)進(jìn)行預(yù)處理，包括離群點(diǎn)剔除和卡爾曼濾波；

2）對(duì)預(yù)處理后的軌跡數(shù)據(jù)進(jìn)行重采樣，使軌跡數(shù)據(jù)的采樣密度基本一致；

3）對(duì)軌跡數(shù)據(jù)點(diǎn)進(jìn)行聚類，得到航路節(jié)點(diǎn)（聚類中心）；

4）對(duì)航路節(jié)點(diǎn)進(jìn)行修正，降低航路節(jié)點(diǎn)的冗余度；

5）根據(jù)軌跡數(shù)據(jù)點(diǎn)的連接關(guān)系，將航路節(jié)點(diǎn)進(jìn)行連接，形成航路幾何數(shù)據(jù)。

1.1 軌跡預(yù)處理

本文涉及的軌跡數(shù)據(jù)由大量離散軌跡數(shù)據(jù)點(diǎn)組成，并存儲(chǔ)在數(shù)據(jù)庫(kù)中。軌跡數(shù)據(jù)點(diǎn)包括軌跡序號(hào)、經(jīng)度、緯度、時(shí)間4 項(xiàng)有效屬性。由于各種因素的影響，在數(shù)據(jù)庫(kù)中存儲(chǔ)的軌跡數(shù)據(jù)存在噪聲和離群點(diǎn)，導(dǎo)致數(shù)據(jù)無(wú)法直接使用。針對(duì)軌跡數(shù)據(jù)中存在的噪聲和離群點(diǎn)，本文提出了如下解決方案。

1.1.1 離群點(diǎn)剔除

采用一種自適應(yīng)閾值法對(duì)離群點(diǎn)進(jìn)行剔除操作。由經(jīng)度、緯度和時(shí)間計(jì)算軌跡點(diǎn)的平均速度，采用基于數(shù)理統(tǒng)計(jì)中的3σ 原則，對(duì)平均速度大于閾值的軌跡點(diǎn)進(jìn)行剔除。

1）閾值初始化。根據(jù)式（1）計(jì)算出整條軌跡各個(gè)采樣點(diǎn)的平均速度。由于中值對(duì)噪聲不敏感，所以本文以平均速度的中值乘以一個(gè)增益系數(shù)η來(lái)初始化閾值。

式中：vi是pi點(diǎn)的速度，E0為初始化閾值，(xi，yi)為軌跡中第i個(gè)軌跡點(diǎn)pi的坐標(biāo)，ti為軌跡點(diǎn)pi對(duì)應(yīng)的時(shí)間。

2）由于當(dāng)前軌跡點(diǎn)的運(yùn)動(dòng)狀態(tài)與相鄰軌跡點(diǎn)的運(yùn)動(dòng)狀態(tài)有一定的相關(guān)性，本文采用滑窗法對(duì)整條軌跡的采樣點(diǎn)進(jìn)行判斷和剔除。該部分分為兩階段：

①窗未滿時(shí)，以初始化閾值為判斷依據(jù)。

式中：pi為軌跡的采樣點(diǎn)，當(dāng)pi=正常點(diǎn)時(shí)，pi進(jìn)入滑窗；否則，將pi從軌跡中刪除。

②窗滿時(shí)，根據(jù)3σ 原則，對(duì)窗內(nèi)軌跡點(diǎn)的平均速度計(jì)算閾值。

式中：E為閾值，meanwin為窗內(nèi)平均速度，varwin為窗內(nèi)速度標(biāo)準(zhǔn)差。

當(dāng)vi>E，pi=離群點(diǎn)，進(jìn)行剔除；否則，窗體滑動(dòng)，將窗體中最舊的點(diǎn)移除窗體，令當(dāng)前點(diǎn)進(jìn)入窗體，直至整條軌跡結(jié)束。

1.1.2 濾噪平滑

采用卡爾曼濾波對(duì)軌跡數(shù)據(jù)進(jìn)行平滑處理，達(dá)到消除軌跡數(shù)據(jù)中部分誤差、平滑軌跡的效果。

1.2 基于軌跡點(diǎn)聚類的航路節(jié)點(diǎn)提取

1.2.1 孤立點(diǎn)剔除

Van Winden 等［15］通過(guò)統(tǒng)計(jì)分析得出結(jié)論：同一道路上的軌跡點(diǎn)到道路中心線的距離近似遵循高斯分布。軌跡預(yù)處理僅降低了軌跡數(shù)據(jù)的噪聲影響，但噪聲仍存在。由于部分軌跡數(shù)據(jù)中會(huì)存在孤立點(diǎn)，偏離航空主干道路的情況，為進(jìn)一步降低對(duì)航路節(jié)點(diǎn)提取的影響，本文采用DBSCAN（Density-Based Spatial Clustering of Applications with Noise）算法［16］對(duì)軌跡數(shù)據(jù)點(diǎn)進(jìn)行深度去噪處理。該算法不僅能夠?qū)壽E數(shù)據(jù)點(diǎn)進(jìn)行聚類，且對(duì)噪聲不敏感，能處理非凸數(shù)據(jù)，較好地識(shí)別孤立點(diǎn)。

定義 孤立點(diǎn)。距離航路中心線較遠(yuǎn)且稀疏的點(diǎn)，如圖1 所示。

圖1 孤立點(diǎn)示意圖Fig.1 Schematic diagram of isolated point

設(shè)軌跡{pi，1，pi，2，…，pi，n}，其中pi，j為第i條軌跡的第j個(gè)軌跡點(diǎn)，軌跡集合I={T1，T2，…，Tm}，軌跡點(diǎn)集合S={p1，1，p1，2，…，pm，n}。

1）對(duì)軌跡點(diǎn)集合S進(jìn)行DBSCAN 聚類并設(shè)置合適的eps和Minpts，可得到聚類結(jié)果S1，S2，…，Sk，S-1，其中，S-1為孤立點(diǎn)集合，其余聚類結(jié)果視為有效點(diǎn)。

2）將孤立點(diǎn)集合S-1中的軌跡點(diǎn)從軌跡集合I中的軌跡T中剔除。

1.2.2 軌跡重采樣

在歷史軌跡數(shù)據(jù)庫(kù)中，軌跡數(shù)據(jù)的采樣時(shí)間間隔不統(tǒng)一，且同一軌跡的采樣點(diǎn)的距離間隔也不同，所以本文采用軌跡重采樣的方法，能消除采樣時(shí)間間隔不統(tǒng)一造成的軌跡點(diǎn)的空間密度不同帶來(lái)的影響，保證各條軌跡數(shù)據(jù)的軌跡點(diǎn)密度基本一致。該步驟的實(shí)現(xiàn)是為了1.2.3 節(jié)中的軌跡點(diǎn)聚類做好數(shù)據(jù)基礎(chǔ)。具體算法如下：

1.2.3 軌跡點(diǎn)聚類

對(duì)數(shù)據(jù)庫(kù)中的歷史數(shù)據(jù)進(jìn)行處理后得到的每條軌跡大致為等距離采樣點(diǎn)，換言之，所有的軌跡數(shù)據(jù)點(diǎn)的密度較均勻，本文采用mini batchK-means 聚類對(duì)所有的軌跡采樣點(diǎn)進(jìn)行聚類，根據(jù)數(shù)據(jù)規(guī)模和覆蓋面積，設(shè)定合適的聚類中心個(gè)數(shù)，經(jīng)過(guò)聚類后，聚類中心能夠較好地反映出軌跡的航路結(jié)構(gòu)，得到航路節(jié)點(diǎn)G={C1，C2，…，Cn}。

1.3 對(duì)航路節(jié)點(diǎn)修正

對(duì)航路節(jié)點(diǎn)進(jìn)行提取后，由于上述聚類方法選取聚類中心個(gè)數(shù)的原因，會(huì)存在航路節(jié)點(diǎn)分裂現(xiàn)象，即兩個(gè)或多個(gè)節(jié)點(diǎn)之間的距離較近，本文采用基于密度的聚類方法，對(duì)航路節(jié)點(diǎn)進(jìn)行合并，得到新的航路節(jié)點(diǎn)集合G′={C′1，C′2，…，C′n′}，并對(duì)軌跡數(shù)據(jù)點(diǎn)pi所屬的節(jié)點(diǎn)簇類進(jìn)行相應(yīng)的修改，如圖2所示。

圖2 航路節(jié)點(diǎn)修正示意圖Fig.2 Schematic diagram of route node amendment

1.4 航路節(jié)點(diǎn)連接

根據(jù)1.3 節(jié)可得到航路節(jié)點(diǎn)，將航路節(jié)點(diǎn)視為聚類中心，那么每條軌跡Ti={pi，1，pi，2，…，pi，n}的軌跡點(diǎn)pi，j(1≤j≤n)都含有各自的所屬類別，本文將軌跡Ti映射到航路節(jié)點(diǎn)集合中，每條軌跡就可由Ti={pi，1，pi，2，…，pi，n}轉(zhuǎn)換為由航路節(jié)點(diǎn)集合中的元素構(gòu)成的序列為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。

軌跡點(diǎn)pi映射到航路節(jié)點(diǎn)集合的原則：

式中：distance(p，C)為軌跡點(diǎn)p和航路節(jié)點(diǎn)C的歐氏距離。

將所有歷史軌跡映射為航路節(jié)點(diǎn)集合中的元素序列后，遍歷數(shù)據(jù)庫(kù)中的所有軌跡，可以得到整個(gè)航路節(jié)點(diǎn)的連接關(guān)系，進(jìn)而得到航路幾何數(shù)據(jù)。

2 實(shí)驗(yàn)與結(jié)果分析

2.1 實(shí)驗(yàn)數(shù)據(jù)集

本文根據(jù)真實(shí)軌跡的運(yùn)動(dòng)特征和數(shù)據(jù)特點(diǎn)，生成仿真數(shù)據(jù)用于實(shí)驗(yàn)。在生成仿真數(shù)據(jù)時(shí)，進(jìn)行添加噪聲、離群點(diǎn)、軌跡斷裂等處理，使仿真數(shù)據(jù)盡可能逼近真實(shí)數(shù)據(jù)。實(shí)驗(yàn)數(shù)據(jù)包括不同噪聲強(qiáng)度的軌跡數(shù)據(jù)，每種噪聲強(qiáng)度下包括4 800條軌跡。實(shí)驗(yàn)數(shù)據(jù)集共有2 296 660 個(gè)軌跡數(shù)據(jù)點(diǎn)，每個(gè)軌跡點(diǎn)包括經(jīng)度、緯度、時(shí)間等屬性。實(shí)驗(yàn)數(shù)據(jù)集如圖3 所示。

圖3 實(shí)驗(yàn)數(shù)據(jù)集Fig.3 Experimental data set

2.2 實(shí)驗(yàn)結(jié)果分析

對(duì)實(shí)驗(yàn)數(shù)據(jù)集進(jìn)行預(yù)處理（如圖4 所示）后，通過(guò)對(duì)孤立點(diǎn)剔除（如圖5 所示），軌跡重采樣、軌跡點(diǎn)聚類和對(duì)航路節(jié)點(diǎn)的修正與連接，得到航路集合數(shù)據(jù)。針對(duì)實(shí)驗(yàn)數(shù)據(jù)集，本文設(shè)定聚類中心個(gè)數(shù)為300，為了對(duì)航路節(jié)點(diǎn)的提取效果進(jìn)行評(píng)價(jià)，本文以仿真生成實(shí)驗(yàn)軌跡數(shù)據(jù)使用的航路結(jié)構(gòu)作為參考，將本文的實(shí)驗(yàn)結(jié)果與參考航路進(jìn)行疊加對(duì)比分析。

圖4 軌跡預(yù)處理前后對(duì)比Fig.4 Comparison of trajectories before and after preprocessing

圖5 孤立點(diǎn)剔除效果Fig.5 Visualization of isolated point removal

2.2.1 定性評(píng)價(jià)

圖6（a）為本文方法的疊加分析圖。從整體結(jié)果來(lái)看，本文方法提取的航路節(jié)點(diǎn)基本覆蓋了參考航路，并且準(zhǔn)確度較高；但對(duì)局部結(jié)果觀察發(fā)現(xiàn)，連接航路節(jié)點(diǎn)的過(guò)程中，部分航路會(huì)出現(xiàn)錯(cuò)誤連接的情況，這主要是由于軌跡數(shù)據(jù)中存在的噪聲的影響，使得軌跡數(shù)據(jù)點(diǎn)向航路節(jié)點(diǎn)集合中映射出現(xiàn)偏差，導(dǎo)致錯(cuò)誤連接。

圖6 兩種方法航路發(fā)現(xiàn)結(jié)果對(duì)比Fig.6 Comparison of results of two route discovery methods

2.2.2 定量評(píng)價(jià)

為了定量評(píng)價(jià)本文方法的實(shí)驗(yàn)結(jié)果，將本文實(shí)驗(yàn)結(jié)果與柵格化方法［8］的實(shí)驗(yàn)結(jié)果分別與參考航路進(jìn)行匹配，然后引用文獻(xiàn)［12］提出的緩沖區(qū)檢測(cè)方法，分別建立10 km、20 km、30 km 的緩沖區(qū)，對(duì)兩種方法提取的航路節(jié)點(diǎn)的覆蓋率α和航路長(zhǎng)度的覆蓋率β進(jìn)行性能評(píng)價(jià)。具體評(píng)價(jià)函數(shù)如下所示：

式中：pi為航路節(jié)點(diǎn)，N為航路節(jié)點(diǎn)總數(shù)，linei為提取的航路中第i段軌跡段在參考航路緩沖區(qū)的部分，LINEi為提取的航路中第i段軌跡段。

如圖7 和圖8 所示，本文方法提取的航路數(shù)據(jù)在節(jié)點(diǎn)覆蓋率和長(zhǎng)度覆蓋率上有較大提高，特別是在10 km 的緩沖區(qū)內(nèi)的實(shí)驗(yàn)結(jié)果。對(duì)單一路段的提取，兩種方法的提取效果比較相近，但對(duì)于兩條航路夾角較小的交匯路口處，柵格化方法易提取成為一條航路；對(duì)“圓”狀航路結(jié)果容易誤提取為點(diǎn)結(jié)構(gòu)。相較于柵格化方法，本文方法能夠在保持提取平滑航路準(zhǔn)確率的基礎(chǔ)上，更好地提取“圓”狀航路結(jié)構(gòu)和夾角較小的兩條航路的航路結(jié)構(gòu)，有較好的抗噪性。

圖7 本文方法與柵格化方法的節(jié)點(diǎn)覆蓋率對(duì)比Fig.7 Comparison of node coverage between two methods

圖8 本文方法與柵格化方法的長(zhǎng)度覆蓋率對(duì)比Fig.8 Comparison of length coverage between two methods

2.3 方法驗(yàn)證

本文方法的有效性已在真實(shí)數(shù)據(jù)集上得到驗(yàn)證，由于真實(shí)數(shù)據(jù)不便公開引用，所以本文對(duì)22 類民航軌跡數(shù)據(jù)進(jìn)行航路提取來(lái)驗(yàn)證本文方法的有效性。軌跡數(shù)據(jù)可視化如圖9（a）所示，航路提取結(jié)果如圖9（b）所示。根據(jù)民航數(shù)據(jù)的航路發(fā)現(xiàn)結(jié)果可知，本文方法能夠較好地發(fā)現(xiàn)航路，但存在一定的局限性，當(dāng)某條航路或航路段上的軌跡數(shù)目稀少時(shí)，該航路或航路段上的軌跡易被檢測(cè)為孤立點(diǎn)進(jìn)行剔除，因此漏掉部分航路。

圖9 民航軌跡數(shù)據(jù)實(shí)驗(yàn)結(jié)果Fig.9 result of civil aviation data set

3 結(jié)語(yǔ)

為了快速準(zhǔn)確地從航空軌跡數(shù)據(jù)中提取航路的幾何數(shù)據(jù)，本文提出了一種基于軌跡點(diǎn)聚類的航路提取方法，利用航空軌跡模擬數(shù)據(jù)進(jìn)行實(shí)驗(yàn)分析，以模擬數(shù)據(jù)的理想航路為參考數(shù)據(jù)，與柵格化方法進(jìn)行定性定量分析評(píng)價(jià)，結(jié)果表明本文方法提取的航路數(shù)據(jù)在節(jié)點(diǎn)覆蓋率和航路長(zhǎng)度覆蓋率上都有明顯提高，而且本文方法易于實(shí)現(xiàn)，有一定的有效性和可靠性。

本文方法提取的航路節(jié)點(diǎn)，可將軌跡數(shù)據(jù)映射為節(jié)點(diǎn)序列，進(jìn)而可將軌跡預(yù)測(cè)問(wèn)題轉(zhuǎn)化為序列到序列的序列預(yù)測(cè)問(wèn)題。為提高預(yù)測(cè)問(wèn)題的正確率，需進(jìn)一步研究和完善：1）軌跡數(shù)據(jù)映射到節(jié)點(diǎn)集合中的更有效的方法；2）各個(gè)航空路段的軌跡數(shù)據(jù)的密度相差較懸殊的情況，本文方法還需要進(jìn)一步改進(jìn)；3）由于噪聲的影響，連接航路節(jié)點(diǎn)時(shí)會(huì)有誤連的情況，下一步需對(duì)連接節(jié)點(diǎn)的方法進(jìn)行進(jìn)一步研究。