基于機(jī)器學(xué)習(xí)的航班靜態(tài)四維軌跡預(yù)測(cè)研究

梁海軍，韓琪聰

(中國(guó)民用航空飛行學(xué)院空中交通管理學(xué)院，四川 廣漢 618307)

1 引言

靜態(tài)四維軌跡預(yù)測(cè)指的是空管部門(mén)在航班飛行前對(duì)航班飛越航路關(guān)鍵點(diǎn)的時(shí)間和高度進(jìn)行估計(jì)?；诤桨囔o態(tài)四維軌跡的估計(jì)結(jié)果，空管部門(mén)可以提前預(yù)知運(yùn)行過(guò)程中空域單元的飛行流量，同時(shí)針對(duì)重大活動(dòng)管制、惡劣天氣等因素對(duì)相關(guān)空域單元的影響，管制部門(mén)也可以提前進(jìn)行航班飛行調(diào)配，以提高空中交通的運(yùn)行效率[1]。同時(shí)，航班運(yùn)行過(guò)程中也可以基于靜態(tài)軌跡估算結(jié)果實(shí)施飛行計(jì)劃與雷達(dá)航跡匹配、飛行沖突和飛行一致性告警等空中交通關(guān)鍵管理手段。因此精確的四維軌跡預(yù)測(cè)是保證飛行安全、維護(hù)空中交通秩序、提高交通效率和經(jīng)濟(jì)性的重要基礎(chǔ)。最初的四維軌跡預(yù)測(cè)方法的基本思想是將整個(gè)航班飛行過(guò)程分為不同階段，按不同階段需要完成的飛行目的建立運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)方程[2]并求解來(lái)完成估算。文獻(xiàn)[3]提出了一種基于飛行剖面劃分的航班軌跡估算。由于這類(lèi)預(yù)測(cè)方法未考慮飛行環(huán)境條件變化對(duì)飛行的影響，且飛行階段的沒(méi)有明確合理劃分，因此預(yù)測(cè)結(jié)果誤差較大。后續(xù)的一些算法引入航空器意圖等因素來(lái)進(jìn)行航班四維軌跡預(yù)測(cè)[4-5]，但該類(lèi)算法的預(yù)測(cè)模型較為簡(jiǎn)單，未能精確的把握航班飛行過(guò)程的隨機(jī)性特點(diǎn)，因此容易造成過(guò)擬合使得預(yù)測(cè)結(jié)果波動(dòng)較大。隨著航班運(yùn)行數(shù)據(jù)的大量存儲(chǔ)[6]，后續(xù)的算法引入了歷史數(shù)據(jù)分析來(lái)進(jìn)行四維軌跡估算[7]，文獻(xiàn)[8]提出了一種基于數(shù)據(jù)融合分析的飛行軌跡估算方法。鑒于上述分析，提出了一種在歷史雷達(dá)或ADS-B(automatic dependency surveillance-broadcast)監(jiān)視數(shù)據(jù)的基礎(chǔ)上，采用概率統(tǒng)計(jì)模型對(duì)飛行過(guò)程中航班的飛行狀態(tài)轉(zhuǎn)移關(guān)系進(jìn)行建模，并通過(guò)歷史數(shù)據(jù)求解最優(yōu)模型參數(shù)，最后以最大后驗(yàn)概率預(yù)測(cè)關(guān)鍵點(diǎn)的飛越時(shí)間和高度。歷史數(shù)據(jù)中包含了航班歷次飛行的航班位置以及環(huán)境因素，是被證明了安全、可行和正確的飛行路線(xiàn)。算法挖掘歷史數(shù)據(jù)規(guī)律而不必對(duì)飛行過(guò)程人為劃分階段進(jìn)行四維軌跡預(yù)測(cè)，結(jié)果更加的精確可信。

2 背景描述

2.1 問(wèn)題描述

現(xiàn)有的空管自動(dòng)化系統(tǒng)以時(shí)間戳10 s對(duì)航空器監(jiān)視數(shù)據(jù)進(jìn)行保存，包含其雷達(dá)監(jiān)視信息、與飛行計(jì)劃匹配信息、預(yù)警告警信息等。在進(jìn)行模型參數(shù)優(yōu)化之前，需要先對(duì)歷史數(shù)據(jù)進(jìn)行處理，將數(shù)據(jù)以航班號(hào)和飛行日期建立搜索索引，且針對(duì)單次飛行的航班數(shù)據(jù)進(jìn)行平滑、去噪和插值等預(yù)處理[9]，并將所有經(jīng)緯度信息轉(zhuǎn)換到直角坐標(biāo)系下。軌跡集采用下式表示：

S={S1，S2，…，ST}

(1)

(2)

其中：式(1)表示按時(shí)間保存的原始估計(jì)數(shù)據(jù)，式(2)表示按航班號(hào)為j飛行日期為i的航班軌跡，共包含t個(gè)位置更新點(diǎn)單個(gè)軌跡點(diǎn)包含x、y、z和t等元素，分別代表直角坐標(biāo)的位置和時(shí)間戳。

靜態(tài)四維軌跡預(yù)測(cè)主要需要解決以下問(wèn)題：

1) 建立航班飛行狀態(tài)轉(zhuǎn)移M，包含未知參數(shù)λ={λ1，λ2，…，λk}；

2) 基于歷史軌跡運(yùn)用合適算法(如EM)迭代優(yōu)化模型未知參數(shù)；

3) 基于最優(yōu)模型預(yù)測(cè)航班靜態(tài)四維軌跡。

2.2 模型描述

基于HMM[10]建立飛行過(guò)程轉(zhuǎn)移模型，HMM通過(guò)描述離散狀態(tài)動(dòng)態(tài)轉(zhuǎn)移和觀(guān)測(cè)值與狀態(tài)之間的關(guān)系來(lái)預(yù)測(cè)對(duì)象未來(lái)時(shí)刻信息，廣泛應(yīng)用于地面交通“停留位置”預(yù)測(cè)，包含如下：

HMM={O，Q，A，B，π}

式中：O表示觀(guān)測(cè)值序列，Q表示隱狀態(tài)序列，A表示隱狀態(tài)轉(zhuǎn)移概率矩陣，B表示觀(guān)測(cè)概率矩陣，π表示隱狀態(tài)初始分布。A中元素aij=p(qt+1=j|qt=i)表示當(dāng)t時(shí)刻狀態(tài)為i時(shí)t+1時(shí)刻狀態(tài)為j的概率，B中元素bij=p(ot=j|qt=i)表示當(dāng)t時(shí)刻隱狀態(tài)為i時(shí)觀(guān)測(cè)值為j的概率。HMM具有如下性質(zhì)(稱(chēng)為馬爾科夫性質(zhì))：

p(qt|qt-1，…，q1，ot-1，…，o1)=p(qt|qt-1)

p(ot|qt，…，q1，ot-1，…，o1)=p(ot|qt)

在建立HMM模型的基礎(chǔ)上，需要解決如下問(wèn)題[12]：

1) 模型參數(shù)學(xué)習(xí)

(3)

上式為EM迭代算法，迭代過(guò)程中未知參數(shù)λ(g+1)>λ(g)總是成立，當(dāng)?shù)鷿M(mǎn)足結(jié)束條件時(shí)的參數(shù)λ(g)即為EM算法下的最優(yōu)參數(shù)。公式中O為歷史監(jiān)視數(shù)據(jù)序列，q∈Q為某一時(shí)刻的航班隱狀態(tài)。

2) 隱狀態(tài)序列估計(jì)

3 飛行狀態(tài)轉(zhuǎn)移建模

采用HMM對(duì)飛行過(guò)程建模的關(guān)鍵是對(duì)模型參數(shù){O，Q，A，B，π}進(jìn)行正確合理描述，既能完整的包含觀(guān)測(cè)值和隱狀態(tài)的所有取值，又能提高算法的效率和預(yù)測(cè)精度。分別對(duì)飛行過(guò)程中位置和高度建模，通過(guò)位置和高度模型預(yù)測(cè)經(jīng)過(guò)關(guān)鍵點(diǎn)的時(shí)間和高度。文章重點(diǎn)介紹位置建模方法，高度建模類(lèi)似。

3.1 觀(guān)測(cè)值建模

歷史軌跡包含航班飛行時(shí)段內(nèi)按周期更新的位置、速度、航向等信息。若以此離散軌跡點(diǎn)對(duì)觀(guān)測(cè)值建模，其數(shù)據(jù)量較大，且相鄰軌跡點(diǎn)屬性由于航空器性能限制不會(huì)產(chǎn)生急劇變化，因此為簡(jiǎn)化模型和提高效率，采用固定網(wǎng)格化地圖對(duì)位置觀(guān)測(cè)值進(jìn)行建模[11]，固定間隔高度對(duì)高度觀(guān)測(cè)值進(jìn)行建模。

從圖1可知，將整個(gè)地圖分為6×8的相同尺寸網(wǎng)格，從左上到右下分別標(biāo)記為(1，1)、(1，2)…(6，8)。劃分網(wǎng)格之后，航班飛行軌跡就可以使用網(wǎng)格標(biāo)記序列表示軌跡。如圖1中軌跡T即可表示為[(1，1)，(2，2)，(3，3)，(3，4)，(3，5)，(4，6)，(5，7)，(5，8)]。顯然，軌跡的觀(guān)測(cè)值序列取決于網(wǎng)格化地圖的劃分方案，當(dāng)劃分不合理時(shí)會(huì)導(dǎo)致軌跡描述缺失問(wèn)題[12]。采用網(wǎng)格化地圖對(duì)觀(guān)測(cè)值建模需要根據(jù)計(jì)算量與預(yù)測(cè)結(jié)果精度進(jìn)行平衡找到最佳參數(shù)。當(dāng)網(wǎng)格尺寸τs(km2)較小時(shí)，建模精度高但計(jì)算量大；反之，在減小計(jì)算量的同時(shí)降低了建模精度，對(duì)應(yīng)的高度層間隔參數(shù)τh(m)對(duì)模型產(chǎn)生相同的影響。建模時(shí)僅需對(duì)飛行計(jì)劃航路包絡(luò)線(xiàn)范圍的地圖及高度層進(jìn)行處理。

圖1 網(wǎng)格化地圖方案示意圖

3.2 隱狀態(tài)建模

本文中HMM隱狀態(tài)才是對(duì)關(guān)鍵點(diǎn)經(jīng)過(guò)時(shí)間和高度最直接的描述，但是從歷史數(shù)據(jù)中不能直接獲得，因此需要通過(guò)與觀(guān)測(cè)值序列建立相應(yīng)關(guān)系進(jìn)行求解。民航航班須在計(jì)劃航路范圍內(nèi)飛行，因此位置隱狀態(tài)集合選擇計(jì)劃航段兩側(cè)偏移10 km的矩形區(qū)域，同時(shí)額外添加停止和偏航2個(gè)狀態(tài)以完整的描述航班位置，圖1中qi即為一個(gè)航段隱狀態(tài)。在高度方向上選擇標(biāo)準(zhǔn)民航飛行高度層作為高度模型隱狀態(tài)。在狀態(tài)轉(zhuǎn)移過(guò)程中，取第一次到達(dá)某一航路段的時(shí)間/高度為預(yù)測(cè)結(jié)果中該航段起始關(guān)鍵點(diǎn)的時(shí)間/高度。

針對(duì)具體飛行計(jì)劃，可以通過(guò)解析AFTN報(bào)文中的航路關(guān)鍵點(diǎn)建立模型隱狀態(tài)集合。圖1中所示的軌跡T包含4個(gè)航路隱狀態(tài)、1個(gè)偏航和1個(gè)停止?fàn)顟B(tài)。

3.3 隱狀態(tài)轉(zhuǎn)移概率矩陣

隱狀態(tài)轉(zhuǎn)移概率矩陣描述了各時(shí)刻航班位置/高度在隱狀態(tài)之間轉(zhuǎn)移的概率。航班飛行狀態(tài)的位置轉(zhuǎn)移的可能性如下：① 當(dāng)前航路段；② 下一航路段；③ 偏航區(qū)域；④ 降落機(jī)場(chǎng)(停止)。狀態(tài)轉(zhuǎn)移之間的關(guān)系是可以描述為“有向圖”，且包含一個(gè)到本狀態(tài)的轉(zhuǎn)移關(guān)系，邊上的權(quán)值為狀態(tài)轉(zhuǎn)移概率。高度方向上狀態(tài)轉(zhuǎn)移包含當(dāng)前以及相鄰的上下高度層。圖2描述了各位置隱狀態(tài)之間轉(zhuǎn)換：

根據(jù)圖2可知，下一時(shí)刻狀態(tài)依然為當(dāng)前狀態(tài)的概率為0.94；為下一狀態(tài)的概率為0.035；為偏航狀態(tài)的概率為0.005；為停止?fàn)顟B(tài)的概率為0.02。每一隱狀態(tài)的轉(zhuǎn)移概率處于矩陣同一行，每行概率和為1表明模型列舉了隱狀態(tài)的所有轉(zhuǎn)移的可能性。

圖2 位置隱狀態(tài)轉(zhuǎn)移有向圖

3.4 觀(guān)測(cè)概率矩陣

觀(guān)測(cè)概率矩陣描述了隱狀態(tài)與觀(guān)測(cè)值之間的概率分布關(guān)系，本文中觀(guān)測(cè)矩陣的意義是對(duì)應(yīng)某一隱狀態(tài)的航班位置/高度分布在某一標(biāo)記的網(wǎng)格化地圖區(qū)域/高度層的概率。同一觀(guān)測(cè)值與隱狀態(tài)的轉(zhuǎn)移概率處于同一行，每行概率和為1表明觀(guān)測(cè)值可能為任一隱狀態(tài)之間存在對(duì)應(yīng)關(guān)系。如圖1所示，假設(shè)t時(shí)刻航班處于隱狀態(tài)q，其觀(guān)測(cè)值概率如表1所示：

表1 狀態(tài)q的觀(guān)測(cè)值概率Table 1 Observation probability of q

3.5 隱狀態(tài)初始分布

根據(jù)馬爾科夫性質(zhì)簡(jiǎn)化模型計(jì)算時(shí)，其中轉(zhuǎn)移概率p(qt|qt-1)和觀(guān)測(cè)概率p(ot|qt)可從模型參數(shù)A和B中獲取，那么還需要在建模時(shí)確定參數(shù)是隱狀態(tài)初始分布p(q1)，本文中即為航班初始位置的分布。由于民航航班在機(jī)場(chǎng)跑道起飛，因此模型中位置隱狀態(tài)初始分布選取以機(jī)場(chǎng)跑道起點(diǎn)為中心的圓形區(qū)域均勻分布；高度隱狀態(tài)初始分布選擇機(jī)場(chǎng)標(biāo)高為均值的均勻分布，分布的待定參數(shù)通過(guò)歷史數(shù)據(jù)進(jìn)行學(xué)習(xí)優(yōu)化。

4 靜態(tài)四維軌跡估算

本章主要針對(duì)位置模型簡(jiǎn)要介紹參數(shù)學(xué)習(xí)和隱狀態(tài)序列預(yù)測(cè)相關(guān)算法[13]，高度模型類(lèi)似。

1) 根據(jù)原始?xì)v史監(jiān)視數(shù)據(jù)信息(航班位置序列和飛行計(jì)劃)完成觀(guān)測(cè)值序列O與隱狀態(tài)序列Q的建模；

2) 基于航班歷史監(jiān)視數(shù)據(jù)軌跡集TS學(xué)習(xí)模型參數(shù)，算法通過(guò)馬爾科夫性質(zhì)簡(jiǎn)化后采用EM算法[7]學(xué)習(xí)模型參數(shù)。

3) 在完成模型參數(shù)學(xué)習(xí)之后，根據(jù)網(wǎng)格化處理之后的觀(guān)測(cè)值序列求出最可能的隱狀態(tài)序列。

取狀態(tài)轉(zhuǎn)移過(guò)程中第一次到達(dá)某個(gè)隱狀態(tài)的時(shí)間/高度作為經(jīng)過(guò)該狀態(tài)起始關(guān)鍵點(diǎn)的時(shí)間和高度；選擇監(jiān)視數(shù)據(jù)更新周期為預(yù)測(cè)周期以保證時(shí)間預(yù)測(cè)精度，當(dāng)更新周期較短時(shí)，可以通過(guò)前向算法來(lái)降低軌跡預(yù)測(cè)算法的時(shí)間復(fù)雜度。令

?i(t)=p(o1，…，ot，qt=i|λ)

(4)

則模型參數(shù)迭代優(yōu)化可簡(jiǎn)化為：

(5)

按照式(4)的定義，有如下遞推關(guān)系：

(6)

其中：a、b為模型參數(shù)A、B中的對(duì)應(yīng)值，k為模型隱狀態(tài)數(shù)量，因此通過(guò)前向算法可將計(jì)算復(fù)雜度從kT降低到k2T，其中T為預(yù)測(cè)周期數(shù)。

5 仿真實(shí)驗(yàn)

本章基于2014年1月1日至2014年12月30日的真實(shí)歷史監(jiān)視數(shù)據(jù)預(yù)測(cè)某一航班12月31日的四維軌跡為例進(jìn)行仿真實(shí)驗(yàn)[14]，首先針對(duì)所有航班數(shù)據(jù)通過(guò)不同的觀(guān)測(cè)值建模參數(shù)τs和τh與預(yù)測(cè)結(jié)果精度之間的關(guān)系確定建模參數(shù)，隨后針對(duì)任一航班預(yù)測(cè)結(jié)果與傳統(tǒng)的運(yùn)動(dòng)學(xué)、動(dòng)力學(xué)和基于歷史數(shù)據(jù)的回歸模型方法的四維軌跡預(yù)測(cè)結(jié)果對(duì)比驗(yàn)證本文算法的有效性。首先引入預(yù)測(cè)結(jié)果誤差評(píng)價(jià)因子[15]：

1) 關(guān)鍵點(diǎn)預(yù)測(cè)誤差

2) 累計(jì)預(yù)測(cè)誤差：

3) 在航班整個(gè)飛行過(guò)程中關(guān)鍵點(diǎn)飛越時(shí)間/高度預(yù)測(cè)誤差的均值和標(biāo)準(zhǔn)差。

5.1 觀(guān)測(cè)值建模參數(shù)

在觀(guān)測(cè)值建模時(shí)，選取不同的網(wǎng)格尺寸以及高度間隔進(jìn)行仿真實(shí)驗(yàn)，比較整個(gè)航班飛行過(guò)程中的經(jīng)過(guò)時(shí)間/高度預(yù)測(cè)均值和方差，以確定最合適的建模參數(shù)τs和τh，在預(yù)測(cè)精度相近的情況下，選擇盡量大的τs和τh以提高算法效率。

根據(jù)實(shí)驗(yàn)結(jié)果可知：

1) 關(guān)鍵點(diǎn)飛越時(shí)間和高度的預(yù)測(cè)誤差均值都隨著網(wǎng)格和高度尺寸變大而變大，但方差總體在一個(gè)較為穩(wěn)定的范圍內(nèi)浮動(dòng)變化；

2) 由表2可知，預(yù)測(cè)時(shí)間誤差均值在網(wǎng)格邊長(zhǎng)小于 4 km時(shí)緩慢增長(zhǎng)，超過(guò)該值后預(yù)測(cè)誤差均值急速增長(zhǎng)，根據(jù)“肘部原則”τs選擇16 km2(邊長(zhǎng)為4 km的正方形)；

3) 由表3可知，預(yù)測(cè)高度均值在高度分層小于30 m時(shí)緩慢增長(zhǎng)，超過(guò)該值后預(yù)測(cè)誤差均值階梯增長(zhǎng)，因此τh選擇30 m。

表2 不同網(wǎng)格化方案下的時(shí)間預(yù)測(cè)均值和標(biāo)準(zhǔn)差Table 2 Means and standard variation of predicted time on different gridded cell

表3 不同高度間隔方案下的高度預(yù)測(cè)均值和標(biāo)準(zhǔn)差Table 3 Means and standard variation of predicted altitude on different altitude separation

5.2 預(yù)測(cè)結(jié)果對(duì)比

選取5.1節(jié)網(wǎng)格尺寸和分層高度間隔參數(shù)進(jìn)行建模、學(xué)習(xí)和預(yù)測(cè)，并與傳統(tǒng)的運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)以及回歸模型預(yù)測(cè)結(jié)果進(jìn)行對(duì)比，其算法預(yù)測(cè)時(shí)間誤差曲線(xiàn)和算法預(yù)測(cè)高度誤差曲線(xiàn)分別如圖3、圖4所示。圖3中縱坐標(biāo)所示的預(yù)測(cè)時(shí)間是自航班起飛時(shí)刻T0開(kāi)始的秒數(shù)。

圖3 算法預(yù)測(cè)時(shí)間誤差曲線(xiàn)

圖4 算法預(yù)測(cè)高度誤差曲線(xiàn)

6 結(jié)論

1) 與傳統(tǒng)的運(yùn)動(dòng)學(xué)、動(dòng)力學(xué)方法和基于歷史數(shù)據(jù)的回歸模型相比：本文提出的方法能更加準(zhǔn)確的把握航班歷史數(shù)據(jù)規(guī)律，不必對(duì)航班飛行過(guò)程劃分階段處理；航班飛越關(guān)鍵點(diǎn)的時(shí)間和高度預(yù)測(cè)結(jié)果也更精確和穩(wěn)定。

2) 以本文提出的方法為基礎(chǔ)，管制部門(mén)實(shí)施的調(diào)控措施能夠提高空中交通運(yùn)行效率。