湯新民 李 騰 陳強超 顧俊偉

(南京航空航天大學民航學院1) 南京 211100) (民航中南地區(qū)空中交通管理局2) 廣州 510000)

0 引 言

在民用飛機的實時運動過程中，外界環(huán)境的擾動將使得航空器偏離預戰(zhàn)術階段生成的能夠有效反應飛行意圖的標稱飛行剖面.因此有必要開發(fā)一種動態(tài)校準4D航跡的算法，不斷追蹤航空器的實時運動狀況，滾動修正航空器未來4D航跡的實時預測值.

NASA的Langley研究中心和Lockfeed-California飛機公司合作，在3D-FMS基礎上進行四維導航研究[1]，著重研究四維航跡的優(yōu)化生成和終端區(qū)的自動化管理技術[2].波音公司也從機載的實時運算方面進行四維航跡算法的研究[3].美國在20世紀90年代開發(fā)的中央塔康自動化系統(tǒng)CTAS是比較突出的一套空中交通管理軟件，它主要針對飛機的進離場次序以及著陸起飛時間開發(fā)的，核心模塊便是四維航跡的生成[4].法國的EEC中心研發(fā)了BADA飛機性能數據庫，而且每年都進行更新，包括200多架不同型號的飛機性能數據，例如燃油流量、飛行包線、飛機質量、發(fā)動機推力、空氣阻力等.BADA為空管系統(tǒng)的航跡仿真和航跡預測計算提供了性能模型[5].

徐琴等[6]建立航空器等速巡航的運動學模型以及等角航跡推測模型，利用ADS-B接收器采集實時航跡數據，運用極大似然法則與牛頓-拉夫遜迭代算法對模型中的巡航地速進行辨識，并利用辨識結果及等角航跡推測模型推算航空器的過點時間.彭瑛等[7]應用大圓航跡和等角航跡原理,融合雷達、電報等動態(tài)數據，研究了動態(tài)航跡推測算法，以及該算法在航空器飛行軌跡預測、航空器過航路點時間預測和空中交通流量動態(tài)統(tǒng)計預測等方面的應用.章濤等[8]提出了一種卡爾曼濾波和拓展卡爾曼濾波聯(lián)合算法辨識運動模型參數的4D航跡預測方法.該方法在等角航跡飛行模型的基礎上，運用KF和EKF聯(lián)合算法辨識航空器的地速，以此計算航空器飛躍未來特征位置的過點時間，余達清[9]用粒子群算法對卡爾曼濾波做出改進，通過優(yōu)化卡爾曼濾波的初始狀態(tài)值和協(xié)方差矩陣提高算法的收斂速度和精度，以彌補初始狀態(tài)值和協(xié)方差矩陣的不確定性對航跡預測造成的不利影響.

上述提到的動態(tài)航跡校準算法大多針對的是等速巡航階段.但航空器的實際飛行是一個相當復雜的過程，多種運動模式共同作用.若用單一的運動模式匹配航空器的實際飛行狀態(tài)，必然造成較大的預測誤差.同時上述算法沒有從駕駛員的實際操縱特性出發(fā)，預測的僅僅是地速，沒有對駕駛員施加給航空器的馬赫數(校準空速)和外界高空大氣的風速、風向深入研究.由此本文提出一種基于多模式交互的短期4D航跡預測算法，用多種運動模式匹配航空器的實際飛行狀態(tài)，對模型中的馬赫數(校準空速)、風速和風向進行辨識.每一種運動模式都對應著一種卡爾曼濾波器，最后的優(yōu)化結果就是各個濾波器輸出的辨識結果的加權平均值[10].

1 構建航空器飛行運動模型

1.1 基于BADA的標稱空速模型

由BADA給出的各個飛行階段的標稱空速剖面為

爬升空速Vclimb:

0≤H≤1 499 ft,Vclimb=CVmin×(Vstall)TO+VdCL,1

1 500 ft≤H≤2 999 ft,Vclimb=CVmin×(Vstall)TO+VdCL,2

3 000 ft≤H≤3 999 ft,Vclimb=CVmin×(Vstall)TO+VdCL,3

4 000 ft≤H≤4 999 ft,Vclimb=CVmin×(Vstall)TO+VdCL,4

5 000 ft≤H≤5 999 ft,Vclimb=CVmin×(Vstall)TO+VdCL,5

6 000 ft≤H≤10 000 ft,Vclimb=min(Vcl,1,250 kn)

10 000 ft≤H≤HTA,Vclimb=Vcl,2

H≥HTA,Vclimb=Macl

(1)

式中：H為航空器的飛行高度；HTA為轉換高度,大約為25 000 ft,見文獻[11]；Vclimb為航空器選定的爬升空速；(Vstall)TO為航空器起飛階段的失速速度；CVmin為修正系數，VdCL,1,VdCL,2,VdCL,3,VdCL,4,VdCL,5為不同飛行高度條件下的爬升空速修正值；Vcl,1,Vcl,2為第一、第二標準爬升校準空速；Macl為標準爬升馬赫數.

巡航空速Vcrulse:

0≤H≤2 999 ft,Vcruise=170 kn

3 000 ft≤H≤5 999 ft,Vcruise=min(Vcr,1,220 kn)

6 000 ft≤H≤13 999 ft,Vcruise=min(Vcr,1,250 kn)

14 000 ft≤H≤HTA,Vcruise=Vcr,2

H≥HTA,Vcruise=Macr

(2)

式中：Vcruise為航空器選定的巡航空速；Macr為標準巡航馬赫數；Vcr,1、Vcr,2為第一、第二標準巡航校準空速.

下降空速Vdescent:

0≤H≤999 ft,Vdescent=CVmin×(Vstall)LD+VdDES,1

1 000 ft≤H≤1 499 ft,Vdescent=CVmin×(Vstall)LD+VdDES,2

1 500 ft≤H≤1 999 ft,Vdescent=CVmin×(Vstall)LD+VdDES,3

2 000 ft≤H≤2 999 ft,Vdescent=CVmin×(Vstall)LD+VdDES,4

3 000 ft≤H≤5 999 ft,Vdescent=min(Vdes,1,220 kn)

6 000 ft≤H≤9 999 ft,Vdescent=min(Vdes,1,250 kn)

6 000 ft≤H≤HTA,Vdescent=Vdes,3

H≥HTA,Vdescent=Mades

(3)

式中：Vdescent為航空器選定的下降空速；Ma為標準下降馬赫數；Vdes,1,Vdes,2為第一、第二標準下降校準空速；Vddes,1,Vddes,2,Vddes,3,Vddes,4為不同飛行高度條件下的下降空速修正值；(Vstall)LD為航空器著陸階段的失速速度.

觀察BADA給出的標稱空速模型，發(fā)現航空器在其提供的標稱高度范圍條件下執(zhí)行的是等校準(CAS)/等馬赫數(Ma)飛行策略.故本文所提出的空速擬合模型基于以下兩種假設條件：

1) 氣象與航向條件保持相對穩(wěn)定 由于現階段主要針對的是短時間內的4D航跡預測，故認為在很短的時間內氣象條件(風速和風向)、航向保持相對穩(wěn)定.

VW(K+1)=VW(K)+σ3×T

θW(K+1)=θW(K)+σ4×T

hd(K+1)=hd(K)+σ5×T

(4)

2) 在標稱高度層內執(zhí)行4D航跡預測 為降低短期4D航跡預測的復雜度，通過ADS-B反饋的飛行高度數據一般都在BADA所提供的標稱高度層內，即不考慮跨越多個標稱高度層的4D航跡預測，暫且認為航空器執(zhí)行等校準空速(CAS)/等馬赫數(Ma)飛行策略.

1.2 狀態(tài)-狀態(tài)轉換方程

參照BADA提出的標準空速模型，又從便于數學處理的原則出發(fā)，認為航空器的運動過程由下述兩種運動模式共同作用產生.

1) 等校準空速(等馬赫數)、垂直方向上的勻速飛行

VH(K+1)=VH(K)+σ1×T

H(K+1)=H(K)+VH(K)×T+σ1×0.5T2

CAS(K+1)=CAS(K)+σ2×

T(M(K+1)=M(K)+σ2×T)

(5)

2) 等校準空速(等馬赫數)、垂直方向上的勻變速飛行

H(K+1)=H(K)+VH(K)×T+

CAS(K+1)=CAS(K)+σ2×

T(M(K+1)=M(K)+σ2×T)

(6)

1.3 狀態(tài)-觀測轉換方程

cos(θW(K+1)-hd(K+1))+β4

TAS(K+1)=39×t(K+1)0.5×0.514×M(K+1)

t(K+1)=T0+ΔT0-1.98(H(K+1)/1 000)

(7)

(8)

1.4 構建觀測轉移矩陣

式(7)～(8)為非線性的轉化關系，需要調用拓展卡爾曼濾波(EKF)來進行線性化處理，生成觀測轉換矩陣B(K+1).拓展卡爾曼濾波(EKF)與卡爾曼濾波(KF)相比，最大的突破之處在于突破了要求線性運動模型的束縛，利用高等數學中的泰勒展開式將非線性的運動模型方程展開，忽略其高階項(二階及二階以上)，只保留一階項和常數項來達到線性化處理的目的，從而將非線性的運動模型納入到濾波處理流程.

(9)

2 基于多模式交互(IMM)的4D航跡濾波與外推

多模式交互算法基本思想是使用多個不同的運動模型分別匹配機動目標的不同運動狀態(tài)，不同模型間的轉移概率是一個馬爾可夫矩陣，上文涉及到非線性的轉化關系，目標狀態(tài)的估計和模型概率的更新使用拓展卡爾曼(EKF)濾波.這樣既克服了使用單模型時一旦目標運動狀態(tài)與模型不符引起的誤差，同時由于各模型之間存在交互，當目標運動模式發(fā)生頻繁跳變時跟蹤精度會大幅度的提高[14].

上文提到航空器的運動過程可以歸結出兩種運動模式，故交互式多模型算法中共有兩個模型,馬爾可夫概率轉移矩陣為P,Pij為從模型i到模型j的轉移概率，且i,j=1,2.

(10)

2.1 輸入交互運算

輸入交互運算是在給定模型上一時刻的狀態(tài)、協(xié)方差估計值并獲得新的量測Z(K+1)之后對模型進行重新初始化運算.首先輸入各子模型的濾波器在K時刻的目標狀態(tài)估計值及其協(xié)方差，即Sj(K)和Pj(K),且i,j=1,2.假設在K時刻匹配的模型是i，在K+1時刻匹配的模型為模型j，計算混合概率得

(11)

[Si(K)-S0j(K)][Si(K)-S0j(K)]′}

(12)

式中：[Si(K)-S0j(K)]′為[Si(K)-S0j(K)]的轉置向量；μi(K)為K時刻匹配模型

2.2 模型概率更新

對于第j個模型，K+1時刻與其匹配的似然函數為

exp{-0.5(dj(K+1))′(Wj(K+1))j-1(dj(K+1))}

(13)

計算各模型的后驗概率得

(14)

2.3 估計融合交互輸出

估計融合是指根據各子模型對應的濾波器輸出的狀態(tài)估計值、協(xié)方差以及經過概率更新后各模型對應的概率，計算出T時刻的目標狀態(tài)的總體估計和總體估計誤差協(xié)方差，分別為

[Sj(K+1)-Sj(K+1)]×

[Sj(K+1)-Sj(K+1)]′}

(15)

2.4 航跡外推

正常情況下，如果沒有獲得ADS-B實時反饋的觀測數據，則需要對航空器展開航跡外推.先調用式(11)～(12)對各運動模型進行歸一化處理，重新初始化各運動模型的狀態(tài)估計值和噪聲協(xié)方差矩陣S0j(m)和P0j(m)(j=1,2).進而利用單模型的狀態(tài)轉換方程(5)～(6)，獲得m+1時刻各個模型的目標狀態(tài)估計值S0j(m+1)和P0j(m+1).唯一與航跡校準階段不同的是模型權重的選擇，具體轉換公式為

(16)

把各模型的預測狀態(tài)值和預測模型概率帶入式(17)便可以得到相互作用多模型的狀態(tài)預測值.

詳細算法流程圖見圖1.

圖1 基于IMM算法的4D航跡預測流程圖

3 案例仿真與討論

根據FlightAware顯示，該航班的機型為A321,BADA給出其在10 000 ft與轉換高度(25 000 ft)之間標稱爬升校準空速為310 kn.調用IMM算法，其辨識后的爬升校準空速穩(wěn)定在300.25 kn.與實際狀況基本符合，存在偏差的原因為本文采用的是拓展卡爾曼濾波(EKF)線性化，忽略了高階項(二次項及二次項以上)，僅僅保留了一階項和常數項，線性化還不夠徹底.

圖2為校準空速、風速、風向的辨識結果圖，由圖2可知，校準空速、風速、風向在第30個觀測點附近趨于收斂，辨識后的風速為15 m/s,辨識后的風向為260°,呈現的是關于上下震蕩且逐漸收斂的圖形.

圖2 辨識結果圖

圖3為垂直速度和地速校準曲線圖，由圖3可知，航空器在爬升過程中垂直(爬升)速度先是逐漸下降，在第34個觀察點附近轉入勻速爬升狀態(tài)，在第58個觀察點附近，爬升速度繼續(xù)衰減.實例表明，勻速爬升和勻變速爬升交替進行，這也說明采用多模式交互算法(IMM)進行濾波的必要性.

圖3 校準曲線圖

通過IMM校準后的垂直速度、地速更加契合觀測值，見表2～3.

表2 算法性能比對表(地速)

表3 算法性能比對表(垂直速度)

由表2～3可知，通過多模式交互算法，與觀測值的誤差平方和、均方誤差明顯降低.但依然與觀測值存在著一定的偏差，原因有二：

1) 選取的觀測點數較多，接近達到81個點，每兩點之間的觀測時間間隔T=5 s.實際的氣象狀況是每一個觀測點的風速和風向都在不斷變化.提出的短時間內的恒定的風速和風向模型與實際略微有點出入.

2) 航空器周圍大氣環(huán)境溫度的計算是基于標準大氣環(huán)境[15]，與實際周圍大氣環(huán)境的溫度必然有所偏差.

利用辨識完成的校準空速、風向、風速大小，針對表1之后的24個航路點的過點時間進行預測，詳細計算方法見文獻[16],預測結果見圖4.本文方法計算出的預計過點時間與實際過點時間的誤差隨著航路點的增多逐漸增大，在前5個航路點內的推測時間誤差可以控制在2 s之內，在前15個航路點內的推測時間誤差可以控制在5 s之內.誤差產生的主要原因是恒定的風速和風向模型過于理想，沒有考慮氣象數據的動態(tài)變化，以及航空器的配載和駕駛員的操縱變化.為了保證預測航跡的準確度，可考慮使用克里金插值算法[17]，用辨識后的風速和風向向下游的航路點進行插值運算，以獲得更準確的氣象信息.

圖4 預計過點時刻和實際過點時刻對比圖

完成IMM算法的歷史航跡校準之后，在缺乏航跡觀測數據的情況下展開航跡外推，推測航空器飛躍下游多個航路點的4D軌跡.勻變速升降模型在歷史航跡校準終端所對應的垂直速度一階導數AVH為-0.302 m/s2，垂直速度為7.06 m/s，飛行高度6 900 m，勻速升降模型在歷史航跡校準終端所對應的垂直速度為7.2 m/s，飛行高度6 930 m.融合多種運動模型之后的總體垂直速度狀態(tài)估計為7.12 m/s，總體飛行速度狀態(tài)估計6 920 m.

針對缺乏歷史觀測數據的11個下游航路點，在不同AVH條件下，生成容納多重運動模式且關于垂直速度和飛行高度的預測剖面圖.

圖5為融合兩種運動模式的總體垂直速度和高度預測剖面圖，AVH為垂直方向勻變速飛行的運動模式在歷史航跡校準的終端所對應的垂直速度加速度(垂直速度的一階導數)，本文所獲得關于AVH仿真結果為-0.302 m/s2，所選取的另外兩值用于對比分析.

圖5 融合兩種運動模式總體預測剖面圖

由圖5可知，隨著AVH在數值上不斷的擴大，垂直速度和飛行高度剖面變化的趨勢趨于變緩.AVH數值上的擴張代表在交互的過程中勻變速升降模式所占據的權重上升.參照BADA性能數據庫提供的校準空速—升降速度的轉換公式，在等校準空速/等馬赫數升降中，航空器的升降速度在爬升段不斷降低，而且隨著飛行高度的增加，升降速度的下降趨于緩慢，這與本文仿真驗證結果基本相符.

4 結束語

本文研究了航空器短期航跡預測問題，利用IMM算法對校準空速(馬赫數)、風向、風速進行了辨識，并利用辨識后的參數提出的全能量模型進行航空器短期4D航跡推測.本文提出的推測方法具備靈活簡便、動態(tài)適應性強的優(yōu)點.下一步的研究重點為將氣象數據精確化，進一步提高航跡預測的準確性，并且建立航空器跨越多個標稱高度層的軌跡預測模型，為優(yōu)化飛行計劃和實施流量管理奠定基礎.