基于交互式多模型的短期4D航跡預(yù)測*

湯新民 李 騰 陳強(qiáng)超 顧俊偉

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

0 引 言

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

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

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

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

1 構(gòu)建航空器飛行運(yùn)動模型

1.1 基于BADA的標(biāo)稱空速模型

由BADA給出的各個飛行階段的標(biāo)稱空速剖面為

爬升空速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為轉(zhuǎn)換高度,大約為25 000 ft,見文獻(xiàn)[11]；Vclimb為航空器選定的爬升空速；(Vstall)TO為航空器起飛階段的失速速度；CVmin為修正系數(shù)，VdCL,1,VdCL,2,VdCL,3,VdCL,4,VdCL,5為不同飛行高度條件下的爬升空速修正值；Vcl,1,Vcl,2為第一、第二標(biāo)準(zhǔn)爬升校準(zhǔn)空速；Macl為標(biāo)準(zhǔn)爬升馬赫數(shù).

巡航空速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為標(biāo)準(zhǔn)巡航馬赫數(shù)；Vcr,1、Vcr,2為第一、第二標(biāo)準(zhǔn)巡航校準(zhǔn)空速.

下降空速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為標(biāo)準(zhǔn)下降馬赫數(shù)；Vdes,1,Vdes,2為第一、第二標(biāo)準(zhǔn)下降校準(zhǔn)空速；Vddes,1,Vddes,2,Vddes,3,Vddes,4為不同飛行高度條件下的下降空速修正值；(Vstall)LD為航空器著陸階段的失速速度.

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

1) 氣象與航向條件保持相對穩(wěn)定 由于現(xiàn)階段主要針對的是短時間內(nèi)的4D航跡預(yù)測，故認(rèn)為在很短的時間內(nèi)氣象條件(風(fēng)速和風(fēng)向)、航向保持相對穩(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) 在標(biāo)稱高度層內(nèi)執(zhí)行4D航跡預(yù)測 為降低短期4D航跡預(yù)測的復(fù)雜度，通過ADS-B反饋的飛行高度數(shù)據(jù)一般都在BADA所提供的標(biāo)稱高度層內(nèi)，即不考慮跨越多個標(biāo)稱高度層的4D航跡預(yù)測，暫且認(rèn)為航空器執(zhí)行等校準(zhǔn)空速(CAS)/等馬赫數(shù)(Ma)飛行策略.

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

參照BADA提出的標(biāo)準(zhǔn)空速模型，又從便于數(shù)學(xué)處理的原則出發(fā)，認(rèn)為航空器的運(yùn)動過程由下述兩種運(yùn)動模式共同作用產(chǎn)生.

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

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) 等校準(zhǔn)空速(等馬赫數(shù))、垂直方向上的勻變速飛行

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)-觀測轉(zhuǎn)換方程

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 構(gòu)建觀測轉(zhuǎn)移矩陣

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

(9)

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

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

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

(10)

2.1 輸入交互運(yùn)算

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

(11)

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

(12)

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

2.2 模型概率更新

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

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

(13)

計算各模型的后驗(yàn)概率得

(14)

2.3 估計融合交互輸出

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

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

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

(15)

2.4 航跡外推

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

(16)

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

詳細(xì)算法流程圖見圖1.

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

3 案例仿真與討論

根據(jù)FlightAware顯示，該航班的機(jī)型為A321,BADA給出其在10 000 ft與轉(zhuǎn)換高度(25 000 ft)之間標(biāo)稱爬升校準(zhǔn)空速為310 kn.調(diào)用IMM算法，其辨識后的爬升校準(zhǔn)空速穩(wěn)定在300.25 kn.與實(shí)際狀況基本符合，存在偏差的原因?yàn)楸疚牟捎玫氖峭卣箍柭鼮V波(EKF)線性化，忽略了高階項(xiàng)(二次項(xiàng)及二次項(xiàng)以上)，僅僅保留了一階項(xiàng)和常數(shù)項(xiàng)，線性化還不夠徹底.

圖2為校準(zhǔn)空速、風(fēng)速、風(fēng)向的辨識結(jié)果圖，由圖2可知，校準(zhǔn)空速、風(fēng)速、風(fēng)向在第30個觀測點(diǎn)附近趨于收斂，辨識后的風(fēng)速為15 m/s,辨識后的風(fēng)向?yàn)?60°,呈現(xiàn)的是關(guān)于上下震蕩且逐漸收斂的圖形.

圖2 辨識結(jié)果圖

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

圖3 校準(zhǔn)曲線圖

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

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

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

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

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

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

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

圖4 預(yù)計過點(diǎn)時刻和實(shí)際過點(diǎn)時刻對比圖

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

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

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

圖5 融合兩種運(yùn)動模式總體預(yù)測剖面圖

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

4 結(jié)束語

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