梁海軍,韓琪聰
(中國民用航空飛行學院空中交通管理學院,四川 廣漢 618307)
靜態(tài)四維軌跡預測指的是空管部門在航班飛行前對航班飛越航路關鍵點的時間和高度進行估計。基于航班靜態(tài)四維軌跡的估計結果,空管部門可以提前預知運行過程中空域單元的飛行流量,同時針對重大活動管制、惡劣天氣等因素對相關空域單元的影響,管制部門也可以提前進行航班飛行調配,以提高空中交通的運行效率[1]。同時,航班運行過程中也可以基于靜態(tài)軌跡估算結果實施飛行計劃與雷達航跡匹配、飛行沖突和飛行一致性告警等空中交通關鍵管理手段。因此精確的四維軌跡預測是保證飛行安全、維護空中交通秩序、提高交通效率和經濟性的重要基礎。最初的四維軌跡預測方法的基本思想是將整個航班飛行過程分為不同階段,按不同階段需要完成的飛行目的建立運動學和動力學方程[2]并求解來完成估算。文獻[3]提出了一種基于飛行剖面劃分的航班軌跡估算。由于這類預測方法未考慮飛行環(huán)境條件變化對飛行的影響,且飛行階段的沒有明確合理劃分,因此預測結果誤差較大。后續(xù)的一些算法引入航空器意圖等因素來進行航班四維軌跡預測[4-5],但該類算法的預測模型較為簡單,未能精確的把握航班飛行過程的隨機性特點,因此容易造成過擬合使得預測結果波動較大。隨著航班運行數據的大量存儲[6],后續(xù)的算法引入了歷史數據分析來進行四維軌跡估算[7],文獻[8]提出了一種基于數據融合分析的飛行軌跡估算方法。鑒于上述分析,提出了一種在歷史雷達或ADS-B(automatic dependency surveillance-broadcast)監(jiān)視數據的基礎上,采用概率統計模型對飛行過程中航班的飛行狀態(tài)轉移關系進行建模,并通過歷史數據求解最優(yōu)模型參數,最后以最大后驗概率預測關鍵點的飛越時間和高度。歷史數據中包含了航班歷次飛行的航班位置以及環(huán)境因素,是被證明了安全、可行和正確的飛行路線。算法挖掘歷史數據規(guī)律而不必對飛行過程人為劃分階段進行四維軌跡預測,結果更加的精確可信。
現有的空管自動化系統以時間戳10 s對航空器監(jiān)視數據進行保存,包含其雷達監(jiān)視信息、與飛行計劃匹配信息、預警告警信息等。在進行模型參數優(yōu)化之前,需要先對歷史數據進行處理,將數據以航班號和飛行日期建立搜索索引,且針對單次飛行的航班數據進行平滑、去噪和插值等預處理[9],并將所有經緯度信息轉換到直角坐標系下。軌跡集采用下式表示:
S={S1,S2,…,ST}
(1)
(2)
其中:式(1)表示按時間保存的原始估計數據,式(2)表示按航班號為j飛行日期為i的航班軌跡,共包含t個位置更新點單個軌跡點包含x、y、z和t等元素,分別代表直角坐標的位置和時間戳。
靜態(tài)四維軌跡預測主要需要解決以下問題:
1) 建立航班飛行狀態(tài)轉移M,包含未知參數λ={λ1,λ2,…,λk};
2) 基于歷史軌跡運用合適算法(如EM)迭代優(yōu)化模型未知參數;
3) 基于最優(yōu)模型預測航班靜態(tài)四維軌跡。
基于HMM[10]建立飛行過程轉移模型,HMM通過描述離散狀態(tài)動態(tài)轉移和觀測值與狀態(tài)之間的關系來預測對象未來時刻信息,廣泛應用于地面交通“停留位置”預測,包含如下:
HMM={O,Q,A,B,π}
式中:O表示觀測值序列,Q表示隱狀態(tài)序列,A表示隱狀態(tài)轉移概率矩陣,B表示觀測概率矩陣,π表示隱狀態(tài)初始分布。A中元素aij=p(qt+1=j|qt=i)表示當t時刻狀態(tài)為i時t+1時刻狀態(tài)為j的概率,B中元素bij=p(ot=j|qt=i)表示當t時刻隱狀態(tài)為i時觀測值為j的概率。HMM具有如下性質(稱為馬爾科夫性質):
p(qt|qt-1,…,q1,ot-1,…,o1)=p(qt|qt-1)
p(ot|qt,…,q1,ot-1,…,o1)=p(ot|qt)
在建立HMM模型的基礎上,需要解決如下問題[12]:
1) 模型參數學習
(3)
上式為EM迭代算法,迭代過程中未知參數λ(g+1)>λ(g)總是成立,當迭代滿足結束條件時的參數λ(g)即為EM算法下的最優(yōu)參數。公式中O為歷史監(jiān)視數據序列,q∈Q為某一時刻的航班隱狀態(tài)。
2) 隱狀態(tài)序列估計

采用HMM對飛行過程建模的關鍵是對模型參數{O,Q,A,B,π}進行正確合理描述,既能完整的包含觀測值和隱狀態(tài)的所有取值,又能提高算法的效率和預測精度。分別對飛行過程中位置和高度建模,通過位置和高度模型預測經過關鍵點的時間和高度。文章重點介紹位置建模方法,高度建模類似。
歷史軌跡包含航班飛行時段內按周期更新的位置、速度、航向等信息。若以此離散軌跡點對觀測值建模,其數據量較大,且相鄰軌跡點屬性由于航空器性能限制不會產生急劇變化,因此為簡化模型和提高效率,采用固定網格化地圖對位置觀測值進行建模[11],固定間隔高度對高度觀測值進行建模。
從圖1可知,將整個地圖分為6×8的相同尺寸網格,從左上到右下分別標記為(1,1)、(1,2)…(6,8)。劃分網格之后,航班飛行軌跡就可以使用網格標記序列表示軌跡。如圖1中軌跡T即可表示為[(1,1),(2,2),(3,3),(3,4),(3,5),(4,6),(5,7),(5,8)]。顯然,軌跡的觀測值序列取決于網格化地圖的劃分方案,當劃分不合理時會導致軌跡描述缺失問題[12]。采用網格化地圖對觀測值建模需要根據計算量與預測結果精度進行平衡找到最佳參數。當網格尺寸τs(km2)較小時,建模精度高但計算量大;反之,在減小計算量的同時降低了建模精度,對應的高度層間隔參數τh(m)對模型產生相同的影響。建模時僅需對飛行計劃航路包絡線范圍的地圖及高度層進行處理。

圖1 網格化地圖方案示意圖
本文中HMM隱狀態(tài)才是對關鍵點經過時間和高度最直接的描述,但是從歷史數據中不能直接獲得,因此需要通過與觀測值序列建立相應關系進行求解。民航航班須在計劃航路范圍內飛行,因此位置隱狀態(tài)集合選擇計劃航段兩側偏移10 km的矩形區(qū)域,同時額外添加停止和偏航2個狀態(tài)以完整的描述航班位置,圖1中qi即為一個航段隱狀態(tài)。在高度方向上選擇標準民航飛行高度層作為高度模型隱狀態(tài)。在狀態(tài)轉移過程中,取第一次到達某一航路段的時間/高度為預測結果中該航段起始關鍵點的時間/高度。
針對具體飛行計劃,可以通過解析AFTN報文中的航路關鍵點建立模型隱狀態(tài)集合。圖1中所示的軌跡T包含4個航路隱狀態(tài)、1個偏航和1個停止狀態(tài)。
隱狀態(tài)轉移概率矩陣描述了各時刻航班位置/高度在隱狀態(tài)之間轉移的概率。航班飛行狀態(tài)的位置轉移的可能性如下:① 當前航路段;② 下一航路段;③ 偏航區(qū)域;④ 降落機場(停止)。狀態(tài)轉移之間的關系是可以描述為“有向圖”,且包含一個到本狀態(tài)的轉移關系,邊上的權值為狀態(tài)轉移概率。高度方向上狀態(tài)轉移包含當前以及相鄰的上下高度層。圖2描述了各位置隱狀態(tài)之間轉換:
根據圖2可知,下一時刻狀態(tài)依然為當前狀態(tài)的概率為0.94;為下一狀態(tài)的概率為0.035;為偏航狀態(tài)的概率為0.005;為停止狀態(tài)的概率為0.02。每一隱狀態(tài)的轉移概率處于矩陣同一行,每行概率和為1表明模型列舉了隱狀態(tài)的所有轉移的可能性。

圖2 位置隱狀態(tài)轉移有向圖
觀測概率矩陣描述了隱狀態(tài)與觀測值之間的概率分布關系,本文中觀測矩陣的意義是對應某一隱狀態(tài)的航班位置/高度分布在某一標記的網格化地圖區(qū)域/高度層的概率。同一觀測值與隱狀態(tài)的轉移概率處于同一行,每行概率和為1表明觀測值可能為任一隱狀態(tài)之間存在對應關系。如圖1所示,假設t時刻航班處于隱狀態(tài)q,其觀測值概率如表1所示:

表1 狀態(tài)q的觀測值概率Table 1 Observation probability of q
根據馬爾科夫性質簡化模型計算時,其中轉移概率p(qt|qt-1)和觀測概率p(ot|qt)可從模型參數A和B中獲取,那么還需要在建模時確定參數是隱狀態(tài)初始分布p(q1),本文中即為航班初始位置的分布。由于民航航班在機場跑道起飛,因此模型中位置隱狀態(tài)初始分布選取以機場跑道起點為中心的圓形區(qū)域均勻分布;高度隱狀態(tài)初始分布選擇機場標高為均值的均勻分布,分布的待定參數通過歷史數據進行學習優(yōu)化。
本章主要針對位置模型簡要介紹參數學習和隱狀態(tài)序列預測相關算法[13],高度模型類似。
1) 根據原始歷史監(jiān)視數據信息(航班位置序列和飛行計劃)完成觀測值序列O與隱狀態(tài)序列Q的建模;
2) 基于航班歷史監(jiān)視數據軌跡集TS學習模型參數,算法通過馬爾科夫性質簡化后采用EM算法[7]學習模型參數。
3) 在完成模型參數學習之后,根據網格化處理之后的觀測值序列求出最可能的隱狀態(tài)序列。
取狀態(tài)轉移過程中第一次到達某個隱狀態(tài)的時間/高度作為經過該狀態(tài)起始關鍵點的時間和高度;選擇監(jiān)視數據更新周期為預測周期以保證時間預測精度,當更新周期較短時,可以通過前向算法來降低軌跡預測算法的時間復雜度。令
?i(t)=p(o1,…,ot,qt=i|λ)
(4)
則模型參數迭代優(yōu)化可簡化為:
(5)
按照式(4)的定義,有如下遞推關系:
(6)
其中:a、b為模型參數A、B中的對應值,k為模型隱狀態(tài)數量,因此通過前向算法可將計算復雜度從kT降低到k2T,其中T為預測周期數。
本章基于2014年1月1日至2014年12月30日的真實歷史監(jiān)視數據預測某一航班12月31日的四維軌跡為例進行仿真實驗[14],首先針對所有航班數據通過不同的觀測值建模參數τs和τh與預測結果精度之間的關系確定建模參數,隨后針對任一航班預測結果與傳統的運動學、動力學和基于歷史數據的回歸模型方法的四維軌跡預測結果對比驗證本文算法的有效性。首先引入預測結果誤差評價因子[15]:
1) 關鍵點預測誤差

2) 累計預測誤差:
3) 在航班整個飛行過程中關鍵點飛越時間/高度預測誤差的均值和標準差。
在觀測值建模時,選取不同的網格尺寸以及高度間隔進行仿真實驗,比較整個航班飛行過程中的經過時間/高度預測均值和方差,以確定最合適的建模參數τs和τh,在預測精度相近的情況下,選擇盡量大的τs和τh以提高算法效率。
根據實驗結果可知:
1) 關鍵點飛越時間和高度的預測誤差均值都隨著網格和高度尺寸變大而變大,但方差總體在一個較為穩(wěn)定的范圍內浮動變化;
2) 由表2可知,預測時間誤差均值在網格邊長小于 4 km時緩慢增長,超過該值后預測誤差均值急速增長,根據“肘部原則”τs選擇16 km2(邊長為4 km的正方形);
3) 由表3可知,預測高度均值在高度分層小于30 m時緩慢增長,超過該值后預測誤差均值階梯增長,因此τh選擇30 m。

表2 不同網格化方案下的時間預測均值和標準差Table 2 Means and standard variation of predicted time on different gridded cell

表3 不同高度間隔方案下的高度預測均值和標準差Table 3 Means and standard variation of predicted altitude on different altitude separation
選取5.1節(jié)網格尺寸和分層高度間隔參數進行建模、學習和預測,并與傳統的運動學和動力學以及回歸模型預測結果進行對比,其算法預測時間誤差曲線和算法預測高度誤差曲線分別如圖3、圖4所示。圖3中縱坐標所示的預測時間是自航班起飛時刻T0開始的秒數。

圖3 算法預測時間誤差曲線

圖4 算法預測高度誤差曲線
1) 與傳統的運動學、動力學方法和基于歷史數據的回歸模型相比:本文提出的方法能更加準確的把握航班歷史數據規(guī)律,不必對航班飛行過程劃分階段處理;航班飛越關鍵點的時間和高度預測結果也更精確和穩(wěn)定。
2) 以本文提出的方法為基礎,管制部門實施的調控措施能夠提高空中交通運行效率。