終端區(qū)航跡簇的中心航跡提取方法研究

趙元棣 ， 孫 禾， 王潔寧， 李 桃

（1. 中國(guó)民航大學(xué)天津市空管運(yùn)行規(guī)劃與安全技術(shù)重點(diǎn)實(shí)驗(yàn)室，天津 300300；2. 中國(guó)民航大學(xué)空中交通管理研究基地，天津 300300；3. 中國(guó)民航大學(xué)圖書館，天津 300300）

終端區(qū)空域系統(tǒng)是整個(gè)空管系統(tǒng)中情況最為復(fù)雜的子系統(tǒng)之一，它是航空器從進(jìn)離場(chǎng)階段到航路飛行階段的過(guò)渡區(qū)域。在終端區(qū)內(nèi)，每個(gè)航班均按照預(yù)先制定的飛行計(jì)劃中的標(biāo)準(zhǔn)進(jìn)離場(chǎng)航線飛行。所謂飛行計(jì)劃，是指航空公司在航班飛行之前，根據(jù)氣象條件、航行情報(bào)、飛機(jī)性能及運(yùn)載重量等情況通過(guò)計(jì)算所制定的標(biāo)準(zhǔn)飛行航線、高度和飛行時(shí)間等。然而，在實(shí)際運(yùn)行中，會(huì)有諸多因素（如：實(shí)時(shí)氣象環(huán)境、終端區(qū)容量、交通管制服務(wù)等）影響航空器的飛行，導(dǎo)致實(shí)際飛行航跡與飛行計(jì)劃所制定的標(biāo)準(zhǔn)航線有所偏差。隨著空中交通流量的不斷增加，終端區(qū)內(nèi)的飛行狀況日益復(fù)雜，給管制員的現(xiàn)場(chǎng)指揮帶來(lái)了巨大的挑戰(zhàn)。因此，研究航班的歷史飛行航跡簇，提取具有代表性的中心航跡并識(shí)別、分析其特征，計(jì)算其與飛行計(jì)劃航跡的偏差，對(duì)于航空公司制定同一航班在未來(lái)的飛行計(jì)劃具有重要的意義，減小了該航班在未來(lái)飛行中出現(xiàn)偏差的可能性，從而有效減輕管制員的工作負(fù)荷。

雷達(dá)監(jiān)視是記錄航空器歷史飛行航跡的一種重要手段，它利用無(wú)線電波發(fā)現(xiàn)目標(biāo)，準(zhǔn)確獲取航空器實(shí)時(shí)的位置、高度、航向等信息。雷達(dá)的掃描周期一般為3～5 s，因此利用雷達(dá)設(shè)備獲取的航空器飛行航跡是一些具有時(shí)序的離散點(diǎn)集；對(duì)于同一航班的歷史航跡簇，從幾何的角度可以看作是一個(gè)在機(jī)場(chǎng)坐標(biāo)系下的三維離散點(diǎn)云集合，可以通過(guò)數(shù)學(xué)的方法對(duì)其進(jìn)行計(jì)算與分析。

計(jì)算幾何是一門利用計(jì)算機(jī)研究數(shù)據(jù)處理和幾何特征的學(xué)科，目前已廣泛應(yīng)用于各種領(lǐng)域。近些年，已有大量學(xué)者成功地將計(jì)算幾何引入民航領(lǐng)域[1-3]；針對(duì)其中的航跡簇聚類與分析問(wèn)題，目前大多都是基于智能計(jì)算的方法[4-6]，而且得到的結(jié)果大多是二維的，在實(shí)際應(yīng)用中具有一定的局限性。骨架提取是計(jì)算幾何中的基礎(chǔ)算法之一，已廣泛應(yīng)用于模型分割[7-8]、目標(biāo)匹配和檢索[9-11]以及曲線曲面重建[12-14]等領(lǐng)域。本文針對(duì)歷史飛行航跡簇，采用骨架提取的方法，提取其中心航跡曲線，并在此基礎(chǔ)上利用平均航向和飛行高度識(shí)別其特征點(diǎn)。首先將歷史飛行航跡簇轉(zhuǎn)化為三維點(diǎn)云模型，接著沿飛行方向進(jìn)行自適應(yīng)采樣得到代表點(diǎn)集合，隨后利用代表點(diǎn)的航向信息建立圖結(jié)構(gòu)并進(jìn)行優(yōu)化得到點(diǎn)云模型的骨架曲線，最后通過(guò)均勻化處理以及光滑處理生成最終的中心航跡曲線。本文通過(guò)對(duì)某終端區(qū)的多條歷史航跡簇真實(shí)數(shù)據(jù)進(jìn)行實(shí)驗(yàn)與分析，驗(yàn)證了該方法的準(zhǔn)確性、有效性和魯棒性。

1 終端區(qū)航跡簇的中心航跡提取方法

1.1 航跡簇生成點(diǎn)云模型

由雷達(dá)設(shè)備接收到的某終端區(qū)飛行航跡數(shù)據(jù)如圖1所示，其中每行數(shù)據(jù)均包括：接收時(shí)間、航班號(hào)、機(jī)場(chǎng)坐標(biāo)系下X坐標(biāo)和Y坐標(biāo)、飛行高度、速度、航向角等信息。由于本文是采用計(jì)算幾何方法對(duì)雷達(dá)數(shù)據(jù)進(jìn)行處理，故只需針對(duì)同一航班在不同日期的飛行航跡提取其幾何相關(guān)信息，即：以機(jī)場(chǎng)歸航雷達(dá)為原點(diǎn)建立的機(jī)場(chǎng)坐標(biāo)系系下的X坐標(biāo)和Y坐標(biāo)，單位為千米(km)；飛行高度，單位為米(m)；以及飛行航向角，單位為度，如圖1中紅框部分所示。

圖1 雷達(dá)監(jiān)視數(shù)據(jù)

按照接收時(shí)間的順序?qū)⑸鲜鰯?shù)據(jù)可視化，可以得到同一航班的飛行航跡簇，如圖2(a)所示，其中每條曲線均代表一條飛行航跡。進(jìn)一步，如果忽略其時(shí)序特征，即可得到飛行航跡簇的三維點(diǎn)云模型P，如圖2(b)所示。需要說(shuō)明的是，由于進(jìn)離場(chǎng)航空器的飛行航跡沿Y軸跨度較大，因此出于美觀的考慮，本文插圖均采用非等比例坐標(biāo)軸。

1.2 沿飛行方向采樣

由于航空器在終端區(qū)內(nèi)最理想的飛行航跡就是飛行計(jì)劃中所制定的標(biāo)準(zhǔn)進(jìn)離場(chǎng)航線，所以由航跡簇生成的三維點(diǎn)云是與標(biāo)準(zhǔn)航線具有相似拓?fù)浣Y(jié)構(gòu)的“管狀”、“無(wú)自交”點(diǎn)云模型；此外，航空器在進(jìn)離場(chǎng)時(shí)速度隨飛行高度降低而減小，導(dǎo)致在雷達(dá)掃描間隔固定的前提下，三維點(diǎn)云模型具有多尺度的特點(diǎn)，即：越靠近X-Y平面的點(diǎn)云密度越大。對(duì)于上述多尺度、無(wú)自交的三維點(diǎn)云 P，本文采用文獻(xiàn)[14]中的方法對(duì)其進(jìn)行自適應(yīng)采樣，具體算法如圖3所示。

圖2 飛行航跡簇及三維點(diǎn)云模型

圖3 自適應(yīng)采樣流程圖

其中Rmin和Rmax為事先給定的閾值。由于航跡簇生成的點(diǎn)云模型在局部大多具有線性結(jié)構(gòu)，因此上述方法使采樣過(guò)程沿著線性伸展區(qū)域，即航空器飛行方向進(jìn)行，從而大大提高了采樣效率。由此方法得到的三維采樣結(jié)點(diǎn)集Q在較好地保持了原始點(diǎn)云模型拓?fù)浣Y(jié)構(gòu)和幾何特征的基礎(chǔ)上，大幅度地減少了三維離散點(diǎn)的數(shù)量，如圖4所示。

圖4 沿飛行方向采樣結(jié)果

1.3 圖的建立與優(yōu)化

利用上述方法得到的每個(gè)采樣結(jié)點(diǎn)qi∈ Q ， 均代表了原始點(diǎn)云模型中的一個(gè)子集Pi? P ，且滿足∪ Pi=P 和 Pi∩Pk= φ，?i≠k。為了準(zhǔn)確表現(xiàn)原始三維點(diǎn)云模型的拓?fù)浣Y(jié)構(gòu)，我們?cè)诓蓸咏Y(jié)點(diǎn)集Q上利用其航向信息建立圖G：對(duì)于兩個(gè)采樣結(jié)點(diǎn)qi,qk∈Q ，當(dāng)它們所對(duì)應(yīng)的原始點(diǎn)云集合Pi和Pk的 K近鄰中含有相同的點(diǎn)，且它們的航向角 hi與hk變化不大（小于給定的閾值）時(shí)，則在它們之間建立一條邊。其中，采樣結(jié)點(diǎn)的航向角設(shè)為其代表點(diǎn)集中點(diǎn)的航向角的平均值。對(duì)于采樣結(jié)點(diǎn)集中的每?jī)蓚€(gè)結(jié)點(diǎn)均采取如上操作，即可得到近似反映原始航跡簇形狀的圖G，如圖5所示。在終端區(qū)中，進(jìn)場(chǎng)航空器在切入五邊進(jìn)近時(shí)的航向角差別較大，因此利用上述方法建立的圖可以避免由于航跡點(diǎn)空間位置相距較近而引起的錯(cuò)誤連接。

圖5 圖G的建立

由圖5可以看出，圖G雖然近似還原出了原始航跡簇的拓?fù)浣Y(jié)構(gòu)，但是在局部含有環(huán)結(jié)構(gòu)，不符合飛行航跡的無(wú)自交特性，因此接下來(lái)利用退化算子對(duì)圖G進(jìn)行優(yōu)化操作：對(duì)于圖G中的環(huán)Ci，計(jì)算其幾何中心作為新的結(jié)點(diǎn)代替該環(huán)，并更新與環(huán)Ci相連的邊集以及其代表點(diǎn)集，如圖6(a)所示。為了在保持圖G的拓?fù)浣Y(jié)構(gòu)的同時(shí)，使圖G的幾何特征更逼近于原始點(diǎn)云模型P，我們從含有度為2的結(jié)點(diǎn)的環(huán)開始應(yīng)用退化算子，直至刪除圖G中所有的環(huán)結(jié)構(gòu)為止。

圖6 圖G的優(yōu)化

在應(yīng)用退化算子進(jìn)行優(yōu)化操作后，圖G已經(jīng)成為一條曲線，如圖6(b～c)所示。然而，在環(huán)結(jié)構(gòu)的退化過(guò)程中，可能會(huì)產(chǎn)生一些多余的短邊，我們利用事先給定的邊長(zhǎng)閾值對(duì)這樣的短邊進(jìn)行識(shí)別并刪除它們，從而得到一條既能準(zhǔn)確反映原始航跡簇的拓?fù)浣Y(jié)構(gòu)，又能近似表示其幾何特征的中心航跡曲線。

1.4 中心航跡的均勻化處理與光滑處理

一般來(lái)講，由上述方法得到的中心航跡曲線的邊長(zhǎng)并不均勻，尤其是在環(huán)結(jié)構(gòu)密集的地方（例如航跡轉(zhuǎn)彎處），我們對(duì)其進(jìn)行均勻化處理，即：在長(zhǎng)邊上插入新的采樣結(jié)點(diǎn)，將其均勻分割成若干短邊；將若干短邊合并成為長(zhǎng)邊，并更新其所代表的點(diǎn)集。此外，由于受航空器性能及空氣動(dòng)力學(xué)限制，航空器在轉(zhuǎn)彎時(shí)的飛行航跡較為光滑，因此我們利用Laplace算子對(duì)中心航跡曲線進(jìn)行光滑處理，使其在幾何上更逼近原始三維點(diǎn)云模型，具體優(yōu)化能量函數(shù)為：

其中，L為圖G所對(duì)應(yīng)的Laplace矩陣，qi′為采樣點(diǎn)qi所代表點(diǎn)集的幾何中心。經(jīng)過(guò)光滑處理后的曲線即為最終的中心航跡曲線，如圖7所示。需要說(shuō)明的是，中心航跡曲線上的點(diǎn)看似分布不均是由非等比例坐標(biāo)軸所致。

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

2.1 中心航跡提取結(jié)果

本文采集到某航班2個(gè)月共60天的進(jìn)場(chǎng)雷達(dá)航跡數(shù)據(jù)，除去非正常進(jìn)近的8條航跡，共有52條正常進(jìn)近的飛行航跡形成航跡簇，應(yīng)用本文方法提取其中心航跡。在實(shí)驗(yàn)中，我們?cè)O(shè)定參數(shù)其中d表示原始點(diǎn)云中每個(gè)點(diǎn)與其鄰居點(diǎn)距離的平均值。

圖8(a～b)所示為由上述數(shù)據(jù)中 25天和 52天的航跡簇生成的三維點(diǎn)云模型，以及提取得到的中心航跡曲線；圖 8(c～d)所示為該終端區(qū)另一航班的飛行航跡簇及與之對(duì)應(yīng)的中心航跡，可以看到本文方法無(wú)論是對(duì)于小規(guī)模的航跡簇（圖8(a)14380個(gè)點(diǎn)，圖8(c)11211個(gè)點(diǎn)），還是對(duì)于大規(guī)模的航跡簇（圖8(b) 27681個(gè)點(diǎn)，圖8(d)27226個(gè)點(diǎn)），均可以有效提取出中心航跡。

需要說(shuō)明的是，中心航跡在轉(zhuǎn)彎點(diǎn)的轉(zhuǎn)彎半徑與航跡簇中不同航跡的轉(zhuǎn)彎半徑分布有關(guān)：轉(zhuǎn)彎半徑分布得越分散，中心航跡的轉(zhuǎn)彎半徑越小。這是由于本文方法在采樣過(guò)程中是將采樣球中點(diǎn)集的重心作為采樣點(diǎn)，而重心的位置依賴于采樣球中點(diǎn)集的分布。因此，本文方法對(duì)于轉(zhuǎn)彎半徑分布密集的航跡簇效果更佳，如圖 8(c～d)。

此外，對(duì)于大規(guī)模的航跡簇，可以通過(guò)增大采樣半徑的方法提高計(jì)算效率。圖9所示為針對(duì)同一個(gè)航跡簇，利用不同采樣半徑進(jìn)行采樣，進(jìn)而提取出的中心航跡曲線。可以看到本文方法對(duì)于采樣半徑不敏感，具有一定的魯棒性。

圖9 利用不同采樣半徑的中心航跡提取

2.2 基于中心航跡的特征點(diǎn)識(shí)別

由本文方法得到的中心航跡曲線中的每一個(gè)采樣點(diǎn)均代表原始點(diǎn)云的一個(gè)子集，因此可進(jìn)一步計(jì)算其平均航向角和平均速度：

其中，ih和vi分別表示原始三維離散點(diǎn)的航向角和速度， c ard ( Pi)表示點(diǎn)集 Pi的元素個(gè)數(shù)。

飛行航跡的特征點(diǎn)包括航向轉(zhuǎn)換點(diǎn)、速度轉(zhuǎn)換點(diǎn)、高度轉(zhuǎn)換點(diǎn)等[15]，它們是4D航跡規(guī)劃與預(yù)測(cè)等多個(gè)問(wèn)題的研究基礎(chǔ)。利用本文方法求得的中心航跡在具有幾何位置信息（X坐標(biāo)、Y坐標(biāo)和飛行高度）的同時(shí)，還可以利用上述方法求得其航向角和飛行速度，從而通過(guò)計(jì)算它們之間的變化，識(shí)別中心航跡上的特征點(diǎn)，圖 10所示為利用中心航跡識(shí)別出的高度轉(zhuǎn)換點(diǎn)（黑色五角星）和航向轉(zhuǎn)換點(diǎn)（粉色五角星）。

圖10 中心航跡的特征點(diǎn)識(shí)別

3 結(jié) 論

本文基于計(jì)算幾何中骨架提取的方法，提出了一種針對(duì)終端區(qū)航跡簇的中心航跡提取方法。在將雷達(dá)航跡數(shù)據(jù)轉(zhuǎn)化為三維點(diǎn)云模型之后，采用沿飛行方向的自適應(yīng)采樣方法對(duì)原始點(diǎn)云進(jìn)行采樣，得到代表點(diǎn)集。充分利用代表點(diǎn)的航向信息建立圖結(jié)構(gòu)，并通過(guò)退化算子進(jìn)行優(yōu)化，形成初始中心航跡曲線。最后先后通過(guò)均勻化處理以及光滑處理，生成最終的中心航跡曲線。在此基礎(chǔ)上，通過(guò)平均航向和平均速度的計(jì)算，識(shí)別中心航跡的特征點(diǎn)。通過(guò)對(duì)真實(shí)航跡簇?cái)?shù)據(jù)進(jìn)行實(shí)驗(yàn)，驗(yàn)證了本文方法的正確性、有效性和魯棒性。本文方法在航跡聚類與分析、4D航跡規(guī)劃與預(yù)測(cè)等領(lǐng)域具有重要的意義和價(jià)值。為了更好地分析終端區(qū)內(nèi)不同航班之間的相互影響，降低發(fā)生飛行沖突的可能性，下一步工作將重點(diǎn)圍繞針對(duì)多航跡簇的中心航跡提取問(wèn)題展開研究。

[1] 騰 鵬,黃 俊,張 斌. 三維航跡的 B 樣條曲線擬合算法[J]. 火力與指揮控制,2007,32(9): 115-118.

[2] 彭建亮,朱 凡,孫秀霞,孫 彪. 基于離散PSO和Voronoi圖的無(wú)人機(jī)航跡規(guī)劃方法研究[C]//第二十六屆中國(guó)控制會(huì)議論文集,湖南: 2007: 804-807.

[3] 李新勝,蘭時(shí)勇,李 綱,肖 朝. 利用 BéZier和 B樣條曲線模擬飛機(jī)航跡的方法[J]. 西南交通大學(xué)學(xué)報(bào),2011,46(6): 1040-1045.

[4] Gariel M,Srivastava A N,Feron E. Trajectory clustering and an application to airspace monitoring [J].IEEE Transactions on Intelligent Transportation Systems,2011,12(4): 1511-1524.

[5] 王增福,潘 泉,郎 林,程詠梅. 基于減法聚類的動(dòng)態(tài)航跡聚類算法[J]. 系統(tǒng)仿真學(xué)報(bào),2009,21(16):5240-5243.

[6] 趙恩來(lái),郝文寧,趙 飛,陳 剛,邵校莎莎. 改進(jìn)的基于密度的航跡聚類算法[J]. 計(jì)算機(jī)工程,2011,37(9): 270-272.

[7] Dey T K,Sun Jian. Defining and computing curve-skeletons with medial geodesic function [C]//Proceedings of the fourth Eurographics symposium on Geometry processing 2006. Aire-la-Ville,Switzerland:Eurographics Association,2006: 143-152.

[8] Au O K-C,Tai C L,Chu H K,Cohen-Or D,Lee T Y.Skeleton extraction by mesh contraction [J]. ACM Transactions on Graphics,2008,27(3): 1-10.

[9] Sundar H,Silver D,Gagvani N,Dickinson S. Skeleton based shape matching and retrieval [C]//Proceedings of the International Conference on Shape Modeling and Applications. Washington DC,USA: IEEE Computer Society,2003: 130-142.

[10] Cornea N D,Demirci M F,Silver D,Shokoufandeh A,Dickinson S J,Kantor P B. 3D object retrieval using many-to-many matching of curve skeletons [C]//Proceedings of the International Conference on Shape Modeling and Applications. Washington DC,USA:IEEE Computer Society,2005: 368-373.

[11] Cao Junjie,Tagliasacchi A,Olson M,Zhang Hao,Su Zhixun. Point cloud skeletons via laplacian-based contraction [C]//Proceedings of the International Conference on Shape Modeling and Applications.Washington DC,USA: IEEE Computer Society,2010:187-197.

[12] Lee I K. Curve reconstruction from unorganized points [J]. Computer Aided Geometric Design,1999,17(2): 161-177.

[13] Tagliasacchi A,Zhang Hao,Cohen Or D. Curve skeleton extraction from incomplete point cloud [J].ACM Transactions on Graphics,2009,28(3): 1-9.

[14] Zhao Yuandi,Cao Junjie,Su Zhixun,Li Zhiyang.Efficient reconstruction of non-simple curves [J].Journal of Zhejiang University-Science C,2011,12(7): 523-532.

[15] 王 超,郭九霞,沈志鵬. 基于基本飛行模型的4D航跡預(yù)測(cè)方法[J]. 西安交通大學(xué)學(xué)報(bào),2009,44(2):296-300.