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

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

1 引言

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

在空中交通流識別方面，學(xué)者們?nèi)〉昧艘欢ǖ难芯砍晒ehm[2]基于二維軌跡數(shù)據(jù)，構(gòu)建基于歐氏距離的相似度矩陣，并利用層次聚類算法對軌跡聚類，識別出關(guān)鍵交通流。Gariel等[3]提出基于K-means和DBSCAN相結(jié)合的交通流識別算法。王超[4]考慮不同類型航空器速度對相似度計(jì)算的影響，對文獻(xiàn)[2]的相似度度量方法進(jìn)行改進(jìn)。王莉莉等人[5]在前人基于航空器位置信息聚類的基礎(chǔ)上，綜合考慮時間、航向和速度等指標(biāo)，構(gòu)建基于LOFC時間窗分割軌跡聚類算法。趙元棣[6]等人結(jié)合重采樣和主成分分析方法，將高維軌跡數(shù)據(jù)映射到低維空間，并基于MeanShift方法建立飛行軌跡聚類方法。文獻(xiàn)[7]采用基于加權(quán)距離和懲罰系數(shù)相結(jié)合的距離度量方法，并提出基于密度的軌跡聚類算法。上述方法盡管能夠?qū)Σ煌沼虻慕煌鬟M(jìn)行識別，均是基于先驗(yàn)知識確定聚類個數(shù)，未完全實(shí)現(xiàn)軌跡聚類的自動化，且聚類速度有待提高。另外，上述方法對于軌跡相似性度量的依賴性較高，容易受到軌跡數(shù)據(jù)規(guī)模等因素的影響，產(chǎn)生錯誤的聚類結(jié)果，可信度較差。

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

2 交通流個數(shù)的自動確定

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

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

1)以機(jī)場基準(zhǔn)點(diǎn)為圓心，構(gòu)建極坐標(biāo)系，并以某半徑r為半徑，確定一個圓形范圍，通過設(shè)定特定的角度將圓形范圍劃分成等大小的多個扇形區(qū)域，如圖2所示。該角度可根據(jù)具體需求來確定，實(shí)現(xiàn)精度可調(diào)節(jié)的目的。

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

3)統(tǒng)計(jì)每個扇形區(qū)域中的軌跡點(diǎn)數(shù)，找到區(qū)間和點(diǎn)跡數(shù)的一一對應(yīng)關(guān)系，繪制曲線圖，如圖2所示。

圖2 極坐標(biāo)構(gòu)建示意圖

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

3 軌跡的相似性度量

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

圖3 軌跡點(diǎn)比對示意圖

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

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

圖4 飛行軌跡重采樣過程

(1)

4 飛行軌跡的譜聚類

目前，比較常用的聚類方法有K-means、DBSCAN和BIRCH方法[8-12]。其中K-means不能解決非凸數(shù)據(jù)；DBSCAN的聚類結(jié)果與參數(shù)設(shè)置有很大的關(guān)系；BIRCH方法的計(jì)算復(fù)雜度較高。譜聚類對數(shù)據(jù)分布的適應(yīng)性更強(qiáng)，聚類效果也很優(yōu)秀，同時聚類的計(jì)算量也小很多，彌補(bǔ)了上述方法的不足。

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

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

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

(2)

其中

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

進(jìn)一步構(gòu)造A的拉普拉斯矩陣

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

(3)

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

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

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

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

圖5 飛行軌跡聚類流程圖

5 實(shí)例分析

以廈門機(jī)場某天256條進(jìn)場飛行軌跡為例，平均每條軌跡包含的軌跡點(diǎn)數(shù)約為3300個，運(yùn)用Matlab軟件編程進(jìn)行聚類計(jì)算。

5.1 交通流個數(shù)的確定

以機(jī)場基準(zhǔn)點(diǎn)為圓心，構(gòu)建極坐標(biāo)系，并以100公里為半徑，確定圓形范圍，對原始軌跡數(shù)據(jù)進(jìn)行裁剪，截取圓形內(nèi)部的軌跡數(shù)據(jù)；將圓形范圍平均分成36個小扇形，每個扇形的角度為10°。

將原始軌跡的經(jīng)緯度坐標(biāo)按墨卡托投影的方法轉(zhuǎn)換為直角坐標(biāo)，再將直角坐標(biāo)轉(zhuǎn)為極坐標(biāo)數(shù)據(jù)，部分極坐標(biāo)數(shù)據(jù)如表1所示。

表1 軌跡點(diǎn)的極坐標(biāo)數(shù)據(jù)

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

5.2 飛行軌跡聚類

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

圖6 基于重采樣數(shù)據(jù)的軌跡聚類結(jié)果

為了驗(yàn)證本文方法的有效性，在沒有對軌跡數(shù)據(jù)重采樣的基礎(chǔ)上，直接對原始軌跡數(shù)據(jù)進(jìn)行聚類，聚類結(jié)果如圖7所示。

圖7 基于原始數(shù)據(jù)的軌跡聚類結(jié)果

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

接下類，本文進(jìn)一步在基于重采樣軌跡數(shù)據(jù)聚類結(jié)果的基礎(chǔ)上，采用對應(yīng)點(diǎn)坐標(biāo)平均取值的方法，計(jì)算各條交通流的代表性軌跡，并在x-z平面繪圖，如圖8所示。

圖8 代表性軌跡提取結(jié)果

6 結(jié)束語

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

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

2)極大值法可有效對交通流的條數(shù)進(jìn)行自動識別，無須認(rèn)為事先設(shè)定，算例結(jié)果表明識別結(jié)果準(zhǔn)確；

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

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