基于軌跡聚類的空中交通流自動識別方法研究

李超群，李善梅，馬維宇，張 程

(中國民航大學空管學院，天津 300300)

1 引言

隨著航空運輸業的快速發展，自改革開放以來，我國民航運輸總周轉量從世界排名第37位躍居世界第2，至2018年，飛機起降架次增長299倍，內地運輸機場數量從78增長到232，旅客吞吐量也是曾經的500倍[1]。空中交通需求的快速增長，使得協同放行、進離場排序、尾隨間隔調配等流量方法得到越來越多的應用。但這些方法僅僅聚焦于實時戰術階段流量調配的研究與應用，主要基于管制經驗，缺乏相關實測運行數據的智能決策分析，未能有效把握空中交通特性變化規律和特征，大大影響流量管理方法的實施效果。航路上航空器密度的逐漸增大，表現出以多航空器聚合形成的交通流特征越發明顯，這就要求基于實際航跡數據，從空中交通流的角度去挖掘空中交通的運行規律。因此，首先需要對空中交通流進行識別。

在空中交通流識別方面，學者們取得了一定的研究成果。Rehm[2]基于二維軌跡數據，構建基于歐氏距離的相似度矩陣，并利用層次聚類算法對軌跡聚類，識別出關鍵交通流。Gariel等[3]提出基于K-means和DBSCAN相結合的交通流識別算法。王超[4]考慮不同類型航空器速度對相似度計算的影響，對文獻[2]的相似度度量方法進行改進。王莉莉等人[5]在前人基于航空器位置信息聚類的基礎上，綜合考慮時間、航向和速度等指標，構建基于LOFC時間窗分割軌跡聚類算法。趙元棣[6]等人結合重采樣和主成分分析方法，將高維軌跡數據映射到低維空間，并基于MeanShift方法建立飛行軌跡聚類方法。文獻[7]采用基于加權距離和懲罰系數相結合的距離度量方法，并提出基于密度的軌跡聚類算法。上述方法盡管能夠對不同空域的交通流進行識別，均是基于先驗知識確定聚類個數，未完全實現軌跡聚類的自動化，且聚類速度有待提高。另外，上述方法對于軌跡相似性度量的依賴性較高，容易受到軌跡數據規模等因素的影響，產生錯誤的聚類結果，可信度較差。

鑒于此，本文提出了空中交通流個數自動確定的極大值方法，并構建基于重采樣和歐氏距離的航空器軌跡相似性度量方法，在此基礎上，提出譜聚類的軌跡聚類算法。本文方法能夠有效克服上述當前航跡聚類方法存在的兩個缺陷，全面實現航跡聚類的自動化，提高聚類準確性的同時，提高聚類的速度。

2 交通流個數的自動確定

空中交通流是大量航空器沿相同航路(段)相繼飛行形成的航班流。由于地面交通流的幾何形狀受道路幾何形狀約束，其幾何形狀基本一致，因此可以直接在某個斷面對交通流的特征信息進行統計。對于空中交通而言，由于航空器飛行受到的空間約束小、受到管制員對飛行軌跡的調控、受對流天氣、導航誤差、飛行技術誤差等不確定因素影響，航空器軌跡的空間分布比較散亂，如圖1所示為廈門機場進近區域某天8:00-24:00所有進場航班軌跡。因此，從散亂的軌跡數據中識別交通流是進行空中交通流特性分析的前提與關鍵。

盡管航空器軌跡的空間分布比較散亂，但沿同一進場程序進場的軌跡具有較高的相似性，軌跡點主要分布在進場程序周圍。鑒于此，本文通過研究軌跡點出現在某一特定區域次數的多少，并將不同區域的次數做成曲線圖，確定其極大值出現的位置與個數，從而自動識別空域交通流的個數，具體計算步驟如下：

1)以機場基準點為圓心，構建極坐標系，并以某半徑r為半徑，確定一個圓形范圍，通過設定特定的角度將圓形范圍劃分成等大小的多個扇形區域，如圖2所示。該角度可根據具體需求來確定，實現精度可調節的目的。

2)由于軌跡數據是以經緯度坐標等記錄航空器的位置信息，不利于進行距離及方位的計算，需要對數據坐標表達方式進行轉換。采用墨卡托投影，將數據從大地經緯度坐標系轉換到直角坐標系，再將數據從直角坐標系轉換到極坐標系。

3)統計每個扇形區域中的軌跡點數，找到區間和點跡數的一一對應關系，繪制曲線圖，如圖2所示。

圖2 極坐標構建示意圖

4)利用Matlab中的[top，location1]=findpeaks(vector)函數求極大值，得到峰值及其個數，從而判斷出空域內交通流的個數。

3 軌跡的相似性度量

目前較常用的相似性度量方法是軌跡點比對法，該方法具有計算簡單，運算速度快等特點。如圖3(a)所示的兩個進場航班的雷達數據的示意圖，每條軌跡上相鄰兩個軌跡點之間的時間差為4秒。沿進場反方向進行兩條軌跡的相應軌跡點之間的比對，逐點計算兩條軌跡的偏差距離，這些偏差距離的平均值便可表征兩條軌跡的相似程度。但考慮到航空器運行的特點，一方面由于不同類型航空器的飛行速度差異，計數相同的軌跡點之間距離在有些情形下并不能準確的描述軌跡之間的差異情況。如圖3(b)中所示的兩條相似的軌跡，飛行速度較快的航空器的軌跡點之間的間距要比速度較慢的航空器軌跡點間距大，這使得相應軌跡點之間的間距實際上成為了軌跡之間的斜距，不能真實反映兩條軌跡的相似程度。另一方面，受氣象，CNS設備，人的不確定性因素的影響，任何兩條軌跡都不可能完全重合，這就導致兩條軌跡所包含的軌跡點數不同，也會影響相似性計算結果的準確性。

圖3 軌跡點比對示意圖

為解決上述兩個問題，本文提出基于重采樣的相似性度量方法，即對劃定區域內的各條軌跡進行重采樣，采樣點數相同。設某條飛行軌跡由n個軌跡點(p1，p2，…，pn)組成。保持首尾軌跡點不變，根據采樣點數(或采樣率)計算參數節點所對應的軌跡點作為重采樣點。例如，采樣點數為m，除保留首尾軌跡點之外，只需計算參數節點為(k-1)/(m-1)(k=2，3，…，m-1)所對應的軌跡點qk即可，則重采樣后的飛行軌跡可表示為(q1，q2，…，qm)，其中，q1=p1，qm=pn。

圖4為重采樣過程的二維示意圖。其中：“○”表示初始軌跡點；“□”表示重采樣點。需要說明的是，該方法可能在重采樣過程中，利用線性插值的方法產生新點，較好地保持了初始軌跡的飛行特征。

圖4 飛行軌跡重采樣過程

(1)

4 飛行軌跡的譜聚類

目前，比較常用的聚類方法有K-means、DBSCAN和BIRCH方法[8-12]。其中K-means不能解決非凸數據；DBSCAN的聚類結果與參數設置有很大的關系；BIRCH方法的計算復雜度較高。譜聚類對數據分布的適應性更強，聚類效果也很優秀，同時聚類的計算量也小很多，彌補了上述方法的不足。

譜聚類是從圖論中演化出來的算法，后來在聚類中得到了廣泛的應用[13-15]。它的主要思想是把所有的數據看做空間中的點，這些點之間可以用邊連接起來。距離較遠的兩個點之間的邊權重值較低，而距離較近的兩個點之間的邊權重值較高，通過對所有數據點組成的圖進行切圖，讓切圖后不同的子圖間邊權重和盡可能的低，而子圖內的邊權重和盡可能的高，從而達到聚類的目的。

本文構建基于譜聚類的飛行軌跡聚類方法，具體步驟如下：

1)基于本文的相似性度量方法，計算兩兩軌跡之間的相似度，構建相似度矩陣A。

(2)

其中

aij=aji，?i=j，aij=1

進一步構造A的拉普拉斯矩陣

L=D-1(D-A)D-1

(3)

2)基于本文的聚類個數確定方法，確定交通流的個數k。

3)提取拉普拉斯矩陣的特征值，并從大到小排列，提取前k個特征值對應的特征向量，構成特征子空間。

4)基于K-means方法，對上述特征子空間進行聚類，聚類個數為k。

本文提出的飛行軌跡聚類方法的主要流程如圖5所示。

圖5 飛行軌跡聚類流程圖

5 實例分析

以廈門機場某天256條進場飛行軌跡為例，平均每條軌跡包含的軌跡點數約為3300個，運用Matlab軟件編程進行聚類計算。

5.1 交通流個數的確定

以機場基準點為圓心，構建極坐標系，并以100公里為半徑，確定圓形范圍，對原始軌跡數據進行裁剪，截取圓形內部的軌跡數據；將圓形范圍平均分成36個小扇形，每個扇形的角度為10°。

將原始軌跡的經緯度坐標按墨卡托投影的方法轉換為直角坐標，再將直角坐標轉為極坐標數據，部分極坐標數據如表1所示。

表1 軌跡點的極坐標數據

接下來統計每個小扇形中的軌跡點數，如表2所示。并計算其極大值，并對低于500的峰值進行舍去，共得到三個峰值，分別為(4574，48274，6703)，分別對應第3號、第24號和第35號扇形，由此可判斷出交通流的個數為3條。

5.2 飛行軌跡聚類

將前述交通流個數的計算結果k=3，輸入到本文軌跡聚類算法中，計算兩兩軌跡之間的相似度，得到相似度矩陣，然后將其轉化為拉普拉斯矩陣，并計算特征向量，最終確定聚類結果，如圖6所示。

圖6 基于重采樣數據的軌跡聚類結果

為了驗證本文方法的有效性，在沒有對軌跡數據重采樣的基礎上，直接對原始軌跡數據進行聚類，聚類結果如圖7所示。

圖7 基于原始數據的軌跡聚類結果

對比分析圖6和圖7，可以看出圖6的聚類效果較好，可看出明顯的3條交通流，而圖7中的軌跡聚類存在誤判現象，未能將三條交通流進行分離，誤判率約為28%。

接下類，本文進一步在基于重采樣軌跡數據聚類結果的基礎上，采用對應點坐標平均取值的方法，計算各條交通流的代表性軌跡，并在x-z平面繪圖，如圖8所示。

圖8 代表性軌跡提取結果

6 結束語

本文從飛行軌跡的特征出發，針對終端區飛行軌跡的聚類問題展開研究，提出一種基于極大值法和重采樣技術的飛行軌跡聚類方法，得到以下結論：

1)基于重采樣的航跡聚類在保留飛行軌跡總體走勢特征的基礎上，對飛行軌跡數據進行降維提高了計算的速度，從而滿足空管自動化對速度的要求；

2)極大值法可有效對交通流的條數進行自動識別，無須認為事先設定，算例結果表明識別結果準確；

3)算例表明，本文方法能夠基于256條飛行軌跡識別出3條交通流，識別效果較好，不存在誤判現象。

4)該方法對于空中交通流運行規律挖掘提供了必要的前提，為后續研究提供準確性保障。