朱代武 劉 豪 黃邦菊 史繼龍 陳澤暉

（1.中國民用航空飛行學(xué)院圖書館 廣漢 618307）（2.中國民用航空飛行學(xué)院空中交通管理學(xué)院 廣漢 618307）

1 引言

隨著航空器數(shù)量的劇增，使得在其運(yùn)行過程中出現(xiàn)空域資源緊張的問題，進(jìn)而造成終端區(qū)和航路擁堵，為此美國和歐洲分別提出下一代空中交通管理系統(tǒng)-NextGen（Next Generation）和 SESAR（Single European Sky ATM Research）［1～3］。上述系統(tǒng)的核心皆是基于航跡的空域運(yùn)行，即航行方式實(shí)現(xiàn)從基于傳感器導(dǎo)航到基于性能導(dǎo)航的轉(zhuǎn)變。同時航跡預(yù)測是輔助管制員決策的關(guān)鍵，因此航跡預(yù)測成為未來空管自動化決策的中堅(jiān)力量。

目前對航跡預(yù)測主要采用以下兩種方式：一是基于卡爾曼波的無參數(shù)估計(jì)。常用的是Lymperopoulos［4］等采用蒙特卡洛法簡介對航跡進(jìn)行評估及Wu［5］等采用的基于歷史數(shù)據(jù)的多元線性回歸法。二是基于航空器性能或氣動學(xué)模型的方法。常見的是Warren等采用航空器運(yùn)動方程，利用混雜系統(tǒng)遞推求解航跡模型。但上述兩種方法經(jīng)常受到航空器機(jī)動性能多樣性、大氣環(huán)境及單一模型濾波的使用限制，因此本文提出基于自適應(yīng)滑動時間窗的四維航跡預(yù)測方法，通過滑動時間窗估計(jì)下一時刻的狀態(tài)矢量，進(jìn)而提高航跡預(yù)測的精度。

2 預(yù)測模型

在航空器實(shí)際運(yùn)行過程中，需通過ADS-B等雷達(dá)設(shè)備及時獲取并進(jìn)行狀態(tài)評估：航班號、位置信息、運(yùn)動狀態(tài)及二次應(yīng)答機(jī)編碼。卡爾曼濾波是一個對動態(tài)序列進(jìn)行線性最小方差估計(jì)的算法，適合用來解決目標(biāo)狀態(tài)跟蹤及位置識別等問題［6～7］。其四維航跡預(yù)測模型的結(jié)構(gòu)示意圖如圖1所示（Visio2015版）。

圖1 四維航跡預(yù)測模型的結(jié)構(gòu)示意圖

圖2 航空器地速分解示意圖

假設(shè)航空器的地速為V，航向?yàn)棣?，對其進(jìn)行速度分解為X軸的速度分量Vx，Y軸的速度分量為Vy即［8］：

卡爾曼濾波根據(jù)系統(tǒng)的狀態(tài)方程和觀測方程，估計(jì)動態(tài)系統(tǒng)下一時刻的狀態(tài)。即：

其中X（m）為系統(tǒng)狀態(tài)矢量，表征在m時刻點(diǎn)的所有的運(yùn)動向量；A（m|m-1）為描述目標(biāo)體運(yùn)動的狀態(tài)轉(zhuǎn)移方程；L（m|m-1）為干擾轉(zhuǎn)移矩陣；w（m）表示運(yùn)動模型的系統(tǒng)噪聲，是p維零均值的白噪聲限；Z（m）為所觀測向量，描述m時刻的系統(tǒng)觀測值；H（m）為所觀測矩陣；v（m）為運(yùn)動估計(jì)過程中所產(chǎn)生的噪聲。

系統(tǒng)噪聲w（m）與觀測噪聲v（m）相互獨(dú)立存在的高斯噪聲，其統(tǒng)計(jì)特性為

其中Q（m）是系統(tǒng)噪聲w（m）的對稱非負(fù)矩陣，R（m）為所觀測到的總噪聲v（m）的對稱方差矩陣。

假設(shè)時間更新狀態(tài)初始分布的均值和協(xié)方差矩陣分別為X0|0=X0和∑0|0=∑0，先驗(yàn)分布p（Xm|H0：m-1）在m時刻的均值和協(xié)方差定義為Xm|m-1和∑m|m-1，后驗(yàn)分布p（Xm|H0：m）在m時刻的均值和協(xié)方差定義為Xm|m和∑m|m［9］，其中先驗(yàn)分布和后驗(yàn)分布又被稱為誤差方差矩陣，可表示為P（m|m-1）。K（m）為濾波增益矩陣，對于m=1，…M，卡爾曼濾波可表示為時間更新和觀測更新的動態(tài)方程［10］：

時間更新方程負(fù)責(zé)推進(jìn)時間軸進(jìn)而預(yù)測下一狀態(tài)的矢量值和誤差協(xié)方差估計(jì)值，為下一刻狀態(tài)構(gòu)建預(yù)評估模型；觀測方程對其矢量及協(xié)方差估計(jì)偏差結(jié)合構(gòu)造校驗(yàn)后的估計(jì)。根據(jù)時間和狀態(tài)更新方程預(yù)測四維航跡從m-1時刻到m時刻的狀態(tài)［10］。

3 滑動時間窗口構(gòu)建

由于卡爾曼濾波本質(zhì)上是將非線性問題線性處理，將會導(dǎo)致濾波單位精度的下降甚至分散。同時超前時間預(yù)測窗體過大會使預(yù)測精度降低，預(yù)測窗口過小航空器無法及時進(jìn)行沖突識別及后續(xù)采取避讓措施，因此根據(jù)飛機(jī)性能、航路環(huán)境建立自適應(yīng)滑動窗體模型，對四維航跡進(jìn)行高精度預(yù)測，其滑動模型如圖3所示。

圖3 自適應(yīng)時間窗體滑動模型

滑動時間窗體可定義為當(dāng)t為超前預(yù)測步數(shù)時，從當(dāng)前時刻m到超前時刻m+t之間的滑動時間窗Tw=［Tm，Tm+t］［11］。通過找尋合適的預(yù)測時間窗體Tw即可求解超前預(yù)測步數(shù)t的過程，就移動目標(biāo)與飛行沖突給出預(yù)測時間窗體的預(yù)測方法。

設(shè)ADS-B傳感器的預(yù)測上限為tmax，超前預(yù)測步數(shù)應(yīng)當(dāng)滿足t≤tmax。由卡爾曼可知加速度隨機(jī)誤差服從正態(tài)分布并且相互獨(dú)立，即εk～N（0，δ2），速度與加速度的估計(jì)值為［vxm，vym，axm，aym］，以y方向的矢量為例，根據(jù)對應(yīng)狀態(tài)運(yùn)動方程可得出：

整理可得：

由正態(tài)分布的可累加性得：

設(shè)uy，m+t為vy，m+t的估計(jì)值，利用［uxm，uym，axm，aym］進(jìn)行航跡矢量信息預(yù)測，可知當(dāng)t＞1時（對X方向同樣適用），并可估算出航空器移動的平均速率u：

令t=tmax，并利用概率分布產(chǎn)生的隨機(jī)矩陣即可得出滑動時間窗體內(nèi)航空器的平均速率。同時采用自適應(yīng)調(diào)整噪聲協(xié)方差P（m|m-1）以減少外界干擾及信號丟失，引入自適應(yīng)矩陣Gm，表述如下：

式中Sm=diag｛s1，s2，s3…sm｝，m為對應(yīng)的航空器矢量維數(shù)，根據(jù)上述系統(tǒng)方程的定義及參數(shù)物理關(guān)系，可得狀態(tài)轉(zhuǎn)移的矩陣A（m|m-1）和觀測矩陣H（m）。

假定航空器的飛行高度滿足Hmin≤h≤hmax，根據(jù)航空器的機(jī)動性能及環(huán)境限制，可構(gòu)造約束條件：

式中Δθ為航向改變角，由于航空器機(jī)動過載限制設(shè)其最大轉(zhuǎn)角為Δθmax；tf為航空器相鄰航跡點(diǎn)的飛行時間；cmin為最短的航跡段長度，即航空器在一個滑動時間窗體內(nèi)的最小直線距離；cmax為最大航程。

4 仿真分析

以ZSNB-ZUUU航路為監(jiān)測對象，以2018年7月份某航班CES9937（A320機(jī)型）的ADS-B實(shí)測數(shù)據(jù)為例，選取寧波櫟社國際機(jī)場的ARP（機(jī)場基準(zhǔn)點(diǎn)）作為空間坐標(biāo)系的原點(diǎn)。部分A320數(shù)據(jù)如表1、表2所示。

表1 CES9937航行時間位置圖

表2 CES9937航行速度高度圖

在仿真實(shí)驗(yàn)中令采樣周期控制在35±5s區(qū)間內(nèi)，觀測噪聲R（n）=0.01×I4×4，Q（n）=I1×1，P（0|0）=8×I6×6，系統(tǒng)噪聲參數(shù)為B（n|n-1）=0.01×I6×1，ADS-B掃描的等時間間隔為2s，因此航空器的航跡由離散的航跡點(diǎn)組成。航班CES9937的離場航跡及巡航航程實(shí)測航跡如圖4所示。

圖4 CES9937航空器ZSNB離場圖

從圖4可知由于信號干擾、阻擋等因素會出現(xiàn)數(shù)據(jù)丟失問題，圓圈部分表示數(shù)據(jù)丟失，按時間間隔2s進(jìn)行m元線性擬合。從圖5可知CES9937航空器的真實(shí)值和觀測值通過卡爾曼波曲線進(jìn)行偏差分析，可知傳統(tǒng)卡爾曼濾波預(yù)測方法對于四維航跡的預(yù)測誤差較大。

圖5 CES9937航空器ZSNB-ZUUU航路圖

采樣周期T=2s，利用自適應(yīng)滑動窗體模型的卡爾曼波并與傳統(tǒng)的卡爾曼波的仿真結(jié)果進(jìn)行仿真，通過利用機(jī)載傳感器1，2對航空器運(yùn)動測量的航跡及濾波預(yù)測結(jié)果進(jìn)行偏差分析，其偏差如圖6所示，同時對自適應(yīng)滑動窗體模型與傳統(tǒng)卡爾曼波X軸和Y軸的預(yù)測誤差進(jìn)行對比［12］。

圖6 基于傳統(tǒng)算法的卡爾曼濾波

由圖6、圖7及圖8可知當(dāng)T=2s時，利用傳統(tǒng)卡爾曼算法進(jìn)行四維航跡預(yù)測時，X軸和Y軸的預(yù)測誤差P（m|m-1）控制在±27m范圍區(qū)間，而利用基于自適應(yīng)滑動時間窗的卡爾曼濾波控制在±22m范圍區(qū)間，對應(yīng)誤差如表3所示。

表3 傳統(tǒng)卡爾曼與基于時間窗卡爾曼對比表

圖7 基于自適應(yīng)滑動時間窗的卡爾曼濾波

圖8 傳統(tǒng)卡爾曼濾波與自適應(yīng)滑動時間窗的卡爾曼濾波對比

按照國際民航組織ICAO對航行通信要求，這里選取±50m作為導(dǎo)航基準(zhǔn)源，區(qū)域概率和航跡的平滑性成正比，區(qū)域概率越大則誤差越大，反之相反。通過對上述數(shù)據(jù)進(jìn)行分析，可知傳統(tǒng)卡爾曼濾波的區(qū)域概率為0.2916，基于滑動時間窗的卡爾曼濾波區(qū)域概率為0.1000，相比可知改進(jìn)后的四維航跡區(qū)域誤差減少65.8%，因此改進(jìn)后的濾波發(fā)散性減弱，航跡的平滑性增強(qiáng)［10～13］，對應(yīng)的濾波增益矩陣K（m）系數(shù)差異減小，按照實(shí)驗(yàn)數(shù)據(jù)對比為表4所示。

表4 傳統(tǒng)卡爾曼與基于時間窗卡爾曼區(qū)域概率對比

以傳統(tǒng)卡爾曼波區(qū)域概率0.2916為全誤差基準(zhǔn)，則改進(jìn)后的誤差概率縮減為原來的（0.2916-0.1000）/0.2916=0.658，對應(yīng)誤差縮減率具體如表5所示。

表5 傳統(tǒng)卡爾曼與基于時間窗卡爾曼區(qū)域誤差對比

由于航空器在實(shí)際飛行過程中，對地航向及速度收到航空器機(jī)動性能的影響不斷變化，預(yù)測誤差往往會上下波動偏置［14］，跟蹤誤差也必然存在波動偏置。通過觀測上述圖中的濾波，改進(jìn)后的滑動時間窗能有效減少εk，在受到外界干擾的情況下其濾波性能夠很快收斂，魯棒性較好［15］。

5 結(jié)語

本文提出基于自適應(yīng)滑動時間窗的四維航跡預(yù)測模型，解決了航空器四維航跡實(shí)時規(guī)劃問題。通過對時間特性賦予的航空器機(jī)動進(jìn)行評估［16］，引入自適應(yīng)矩陣Gm保持卡爾曼濾波中的混合高斯模型適應(yīng)外界干擾的多模態(tài)特點(diǎn)［17］，并在此基礎(chǔ)上利用矢量因素及設(shè)定的誤差矩陣P（m|m-1）求解實(shí)時四維航跡［18～20］。仿真結(jié)果證明基于自適應(yīng)滑動時間窗的四維航跡預(yù)測模型是有效的，利用機(jī)載傳感器及ADS-B數(shù)據(jù)能夠?qū)嵤┚珳?zhǔn)四維航跡預(yù)測。