基于改進(jìn)滑動(dòng)窗多項(xiàng)式擬合法的航空器4D航跡預(yù)測

張建祥，甘旭升，周志靖， 楊捷

(1.西京學(xué)院 理學(xué)院，西安 710123) (2.空軍工程大學(xué) 空管領(lǐng)航學(xué)院，西安 710051)

0 引 言

航跡是指航空器在飛行過程中按時(shí)間順序所經(jīng)歷的全部或部分空間位置點(diǎn)的集合，將時(shí)間與各個(gè)空間點(diǎn)位一一對應(yīng)后的組合被稱為4D航跡。航跡預(yù)測是指根據(jù)航空器初始飛行計(jì)劃、歷史航跡信息與經(jīng)驗(yàn)信息，采用相應(yīng)預(yù)測方法，對航空器未來時(shí)間段內(nèi)將要產(chǎn)生的航跡點(diǎn)進(jìn)行預(yù)測。高精度的4D航跡預(yù)測是無人機(jī)飛行防相撞空中威脅態(tài)勢預(yù)警技術(shù)的基礎(chǔ)，是提高空中交通管理效率與安全水平的重要技術(shù)手段。

卡爾曼濾波是由 R.E.Kalman[1]于1960年推導(dǎo)出來的一種最小方差最優(yōu)估計(jì)算法，主要用于導(dǎo)航系統(tǒng)多傳感器信息融合處理和目標(biāo)狀態(tài)估計(jì)。1998年，L.Tong等[2]首次應(yīng)用卡爾曼濾波算法解決航跡預(yù)測問題，對相鄰時(shí)刻目標(biāo)位置及其瞬時(shí)速度進(jìn)行預(yù)測，從一定程度上改善了預(yù)測精度，但在實(shí)際應(yīng)用中存在非線性系統(tǒng)的適用性問題，預(yù)測性能還有較大提升空間。之后，又有諸多研究人員先后提出改進(jìn)型的基于殘差均值的交互式多模型跟蹤濾波算法處理航跡預(yù)測問題，該濾波算法對機(jī)動(dòng)目標(biāo)追蹤預(yù)測有較高精度，但需要做大量的數(shù)學(xué)建模工作，例如航空器運(yùn)動(dòng)模型、運(yùn)行環(huán)境模型等，并最終融合各種模型，形成航跡預(yù)測模型，工作量非常巨大[3-6]。此外，引入數(shù)據(jù)挖掘方法對航空器歷史航跡數(shù)據(jù)進(jìn)行處理、分析，挖掘數(shù)據(jù)中隱含的模式，從而預(yù)測航空器未來位置，并已在航空器4D航跡預(yù)測領(lǐng)域得到應(yīng)用，比較成熟的方法，例如神經(jīng)網(wǎng)絡(luò)算法，存在局部收斂問題和網(wǎng)絡(luò)結(jié)構(gòu)確定問題。美國艾姆斯研究中心與美國管制員協(xié)會(huì)聯(lián)合研發(fā)了中央終端區(qū)管制自動(dòng)化系統(tǒng)[7]，其航跡預(yù)估模塊是核心部分，可生成空域內(nèi)所有航空器未來30 min內(nèi)的航跡。PATs[8]是歐洲開發(fā)的航跡預(yù)測工具，可預(yù)先評估所發(fā)出的管制指令對航空器的影響，提高歐洲地區(qū)空中交通管理效率與安全水平。國內(nèi)的航空器4D航跡預(yù)測研究起步較晚。張軍峰等[9]提出了基于連續(xù)動(dòng)態(tài)模型與離散動(dòng)態(tài)模型的4D航跡預(yù)測方法，并將其用于離場航跡預(yù)測，但該方法復(fù)雜度較高，無法滿足實(shí)時(shí)性需求；王建忠等[10]根據(jù)點(diǎn)融合進(jìn)近的特點(diǎn)，提出航空器4D航跡預(yù)測方法，并采用遺傳算法進(jìn)行求解，但該方法僅能處理簡單問題，無法預(yù)測多跑道航跡的復(fù)雜情況；張振興等[11]提出基于實(shí)時(shí)反饋長短期記憶神經(jīng)網(wǎng)絡(luò)的4D航跡預(yù)測模型，但該方法確定網(wǎng)絡(luò)參數(shù)存在較大難度，所建模型通常無法保證精度。

本文提出采用滑動(dòng)窗多項(xiàng)式擬合法解決航空器4D航跡預(yù)測問題。對滑動(dòng)窗多項(xiàng)式擬合法進(jìn)行改進(jìn)，構(gòu)建4D航跡預(yù)測模型，力求給出一種精度高、簡單易行且計(jì)算簡便的航跡預(yù)測方法，并通過實(shí)例驗(yàn)證其有效性和可行性。

1 滑動(dòng)窗多項(xiàng)式擬合預(yù)測原理

曲線擬合的具體作法：針對給定數(shù)據(jù)點(diǎn){Xi|i=0,1,…,n}，在給定的函數(shù)類Φ中，求f(ti)∈Φ，使誤差ri=f(ti)-Xi(i=0,1,…,n)的平方和最小，即

(1)

從幾何意義上講，就是找尋出一條曲線y=f(t)，使之與給定點(diǎn){Xi|i=0,1,…,n}的距離平方和為最小。函數(shù)f(t)稱為擬合函數(shù)或最小二乘解，求擬合函數(shù)f(t)的方法稱為曲線擬合的最小二乘法。

(2)

A=(PTP)-1PTX

(3)

易證式(3)中(PTP)為對稱正定矩陣，故存在唯一解。求解式(3)可得擬合多項(xiàng)式

(4)

由式(4)可對(tn+dt)時(shí)刻的目標(biāo)航跡進(jìn)行預(yù)測，即采用多項(xiàng)式擬合法預(yù)測航跡：

(5)

固定式(5)中的n值，即預(yù)測模型取固定數(shù)量的數(shù)據(jù)點(diǎn)，并隨著周期更新，不斷接收數(shù)據(jù)點(diǎn)，剔除相同數(shù)量的舊數(shù)據(jù)點(diǎn)，即在多項(xiàng)式擬合法中應(yīng)用滑動(dòng)窗思想，模型將具有動(dòng)態(tài)預(yù)測能力。

2 滑動(dòng)窗多項(xiàng)式擬合法的改進(jìn)

為了最大程度上利用歷史數(shù)據(jù)信息，提高算法預(yù)測性能，對傳統(tǒng)滑動(dòng)窗多項(xiàng)式擬合法提出如下改進(jìn)措施。

2.1 構(gòu)建合適多項(xiàng)式擬合方程

傳統(tǒng)滑動(dòng)窗多項(xiàng)式擬合法的優(yōu)點(diǎn)在于采用遞推法最大程度地利用有限的歷史數(shù)據(jù)信息，而本文對滑動(dòng)窗多項(xiàng)式擬合法的改進(jìn)之處則是進(jìn)一步強(qiáng)化了該優(yōu)點(diǎn)。傳統(tǒng)滑動(dòng)窗多項(xiàng)式擬合法基于歷史數(shù)據(jù)組構(gòu)建了一個(gè)多項(xiàng)式擬合方程，依據(jù)該多項(xiàng)式擬合方程計(jì)算所有預(yù)測值，但每個(gè)預(yù)測值與歷史數(shù)據(jù)組的時(shí)間距離均不同，僅僅依靠同一個(gè)多項(xiàng)式擬合方程，并不能完全表達(dá)出每個(gè)未來值受歷史數(shù)據(jù)組的影響程度。本文在同時(shí)預(yù)測多個(gè)連續(xù)未來值時(shí)，為每個(gè)預(yù)測值構(gòu)造了更為合適的多項(xiàng)式擬合方程，以便更準(zhǔn)確地體現(xiàn)歷史數(shù)據(jù)組對每個(gè)未來值的不同影響，將得到比傳統(tǒng)滑動(dòng)窗多項(xiàng)式擬合法更高的預(yù)測精度。

為每一個(gè)預(yù)測值構(gòu)建更合適的多項(xiàng)式擬合方程的關(guān)鍵在于挑選合適的歷史數(shù)據(jù)組，本文基于每個(gè)預(yù)測值與當(dāng)前值的時(shí)間距離來挑選歷史數(shù)據(jù)組，充分考慮多項(xiàng)式擬合法的時(shí)間序列性質(zhì)。在實(shí)際應(yīng)用中，假設(shè)當(dāng)前值為Xk，滑動(dòng)窗寬度為n，即預(yù)測模型中的歷史數(shù)據(jù)個(gè)數(shù)為n，預(yù)測當(dāng)前值之后1,2,…,q時(shí)刻的未來值，針對當(dāng)前值之后第s(1≤s≤q)時(shí)刻的未來值進(jìn)行預(yù)測時(shí)，其與當(dāng)前值的時(shí)間距離為s，則其歷史數(shù)據(jù)組的選擇應(yīng)為{Xk-s(n-1),Xk-s(n-2),…,Xk-s,Xk}，再依據(jù)式(5)，即可為每個(gè)預(yù)測值提供更合適的多項(xiàng)式擬合方程。

2.2 滑動(dòng)窗多項(xiàng)式參數(shù)自適應(yīng)

本文所提出的滑動(dòng)窗多項(xiàng)式參數(shù)自適應(yīng)主要是指擬合多項(xiàng)式階數(shù)與滑動(dòng)窗長度的自適應(yīng)。

2.2.1 擬合多項(xiàng)式的階數(shù)自適應(yīng)

根據(jù)多項(xiàng)式的最小二乘擬合原理可知，在計(jì)算預(yù)測模型時(shí)，首先要確定擬合多項(xiàng)式的階數(shù)m。要獲得高精度的航跡預(yù)測值，必須判明目標(biāo)的運(yùn)動(dòng)狀態(tài)，使擬合多項(xiàng)式的階數(shù)與其相適應(yīng)。本文將目標(biāo)運(yùn)動(dòng)模式劃分為類直線運(yùn)動(dòng)與曲線運(yùn)動(dòng)，設(shè)定目標(biāo)處于類直線運(yùn)動(dòng)模式與曲線運(yùn)動(dòng)模式時(shí)，對應(yīng)的擬合多項(xiàng)式階數(shù)m分別為1和2。

對于三維空間中目標(biāo)的運(yùn)動(dòng)模式是否存在類直線與曲線之間的突變以及如何進(jìn)行突變的判定問題，將其簡化到三個(gè)單維空間之中進(jìn)行解決，根據(jù)目標(biāo)在單個(gè)維度之中的航跡點(diǎn)繪制時(shí)間-位置折線圖，取相鄰折線所形成角度的大小以及角度的連續(xù)變化趨勢為判定依據(jù)。

在目標(biāo)航跡序列中取X軸上相鄰的四個(gè)航跡點(diǎn)為xk-3,xk-2,xk-1,xk，由于具有相同的時(shí)間間隔，根據(jù)上述四個(gè)航跡點(diǎn)可直接繪制成時(shí)間-位置折線圖，如圖1所示，α1為目標(biāo)在X軸上航跡點(diǎn)xk-3,xk-2,xk-1形成相鄰折線的夾角，α2為目標(biāo)在X軸上航跡點(diǎn)xk-2,xk-1,xk形成相鄰折線的夾角。

根據(jù)向量夾角公式可得出圖1中兩個(gè)連續(xù)夾角的計(jì)算方程為[13]

(6)

式中：arccos為三角函數(shù)中反余弦符號；‖‖2為向量的2-范數(shù)符號[14]。

圖1 時(shí)間-位置折線圖Fig.1 Time-position line chart

可通過設(shè)定一個(gè)參數(shù)λ來進(jìn)行目標(biāo)運(yùn)動(dòng)模式與相關(guān)突變情況的判定：

若目標(biāo)航跡滿足α1>λ,α2>λ，則目標(biāo)在X軸上處于曲線運(yùn)動(dòng)模式，擬合多項(xiàng)式階數(shù)取m=2；若目標(biāo)航跡滿足α1>λ,α2≤λ，則目標(biāo)在X軸上是從曲線運(yùn)動(dòng)模式轉(zhuǎn)為類直線運(yùn)動(dòng)模式，擬合多項(xiàng)式階數(shù)取m=1；若目標(biāo)航跡滿足α1≤λ,α2>λ，則目標(biāo)在X軸上是從類直線運(yùn)動(dòng)模式轉(zhuǎn)為曲線運(yùn)動(dòng)模式，擬合多項(xiàng)式階數(shù)取m=2；若目標(biāo)航跡滿足α1<λ,α2<λ，則目標(biāo)在X軸上是處于類直線運(yùn)動(dòng)模式，擬合多項(xiàng)式階數(shù)取m=1。

目標(biāo)在Y軸與Z軸上的運(yùn)動(dòng)模式與相關(guān)突變情況的判定，與目標(biāo)在X軸上的判定原理相同，此處不再贅述。

2.2.2 多項(xiàng)式滑動(dòng)窗長度自適應(yīng)

因?yàn)轭A(yù)測是基于有限的歷史數(shù)據(jù)進(jìn)行的，所以必須選擇所使用時(shí)間窗寬度。如果時(shí)間窗寬度較大，即歷史數(shù)據(jù)太多，則預(yù)測算法的運(yùn)算量加大，導(dǎo)致預(yù)測緩慢，且如果目標(biāo)運(yùn)動(dòng)軌跡變化較快，將導(dǎo)致實(shí)際目標(biāo)運(yùn)動(dòng)軌跡與時(shí)間多項(xiàng)式函數(shù)之間的擬合效果惡化，則預(yù)測偏差增加；如果時(shí)間窗寬度太小，則數(shù)據(jù)量過少，所計(jì)算出的多項(xiàng)式函數(shù)將不能描述目標(biāo)的運(yùn)動(dòng)特性，預(yù)測效果也不理想。因此，確定滑動(dòng)窗寬度時(shí)，一方面要考慮不能使預(yù)測模型計(jì)算太多的歷史數(shù)據(jù)，即滑動(dòng)窗寬度不可太寬，另一方面則要兼顧目標(biāo)運(yùn)動(dòng)軌跡變化的快慢，即須考慮所選取的擬合多項(xiàng)式階數(shù)。

綜上所述，為了進(jìn)一步減少擬合誤差，提高預(yù)測精度，依據(jù)多項(xiàng)式擬合原理得出擬合多項(xiàng)式滑動(dòng)窗長度n與階數(shù)m之間的關(guān)系式(式(7))，而且在低階數(shù)多項(xiàng)式擬合法的位置預(yù)測研究中，李亞寧[12]已通過仿真證明了其合理性。

n=m+1

(7)

3 多項(xiàng)式擬合法航跡預(yù)測實(shí)現(xiàn)

采用改進(jìn)的滑動(dòng)窗多項(xiàng)式擬合法對航空器航跡進(jìn)行預(yù)測，其過程是先在直角坐標(biāo)系中X,Y,Z軸上分別對航空器航跡進(jìn)行預(yù)測，而后形成對航空器三維空間航跡的總體預(yù)測。具體步驟如下：

(1) 按照改進(jìn)方法在X,Y,Z軸上為每個(gè)預(yù)測值挑選合適的歷史數(shù)據(jù)組。

(2) 依據(jù)式(6)及目標(biāo)運(yùn)動(dòng)模式與相關(guān)突變情況的判定結(jié)果，確定X,Y,Z軸上的擬合多項(xiàng)式階數(shù)mx,my,mz，滑動(dòng)窗長度nx,ny,nz。航空器在X,Y,Z軸上的運(yùn)動(dòng)軌跡多項(xiàng)式擬合函數(shù)為

(8)

(9)

式中：T為更新周期。

求得式(9)中X,Y,Z軸上航空器航跡預(yù)測值，即可得航空器三維空間中航跡總體預(yù)測值。

(4) 采用遞推算法，重新為每個(gè)預(yù)測值挑選合適的歷史數(shù)據(jù)組，并利用滑動(dòng)窗口不斷更新最新的X,Y,Z軸上航空器航跡數(shù)據(jù)，重新確定X,Y,Z軸上的擬合多項(xiàng)式階數(shù)mx,my,mz與滑動(dòng)窗寬度nx,ny,nz，求得基于線性最小方差估計(jì)的多項(xiàng)式預(yù)測模型，預(yù)測下一周期航空器的航跡。

4 實(shí)驗(yàn)仿真與分析

在MATLAB 2014a計(jì)算環(huán)境中，分析改進(jìn)滑動(dòng)窗多項(xiàng)式擬合法預(yù)測航跡的性能，并與傳統(tǒng)滑動(dòng)窗多項(xiàng)式擬合法進(jìn)行對比。由于后者無法自適應(yīng)調(diào)整多項(xiàng)式階數(shù)與滑動(dòng)窗長度，在仿真分析中，將根據(jù)階數(shù)與滑動(dòng)窗長度不同，設(shè)置多個(gè)對照組。

4.1 數(shù)據(jù)來源

實(shí)驗(yàn)使用數(shù)據(jù)來自于通過六階貝塞爾曲線法構(gòu)建的入侵機(jī)航跡數(shù)據(jù)，即入侵機(jī)相對無人機(jī)的航跡數(shù)據(jù)，如圖2～圖5所示。在所有的原始數(shù)據(jù)中，將其中0～99 s的數(shù)據(jù)只作為歷史數(shù)據(jù)使用，預(yù)測航跡為100～196 s中每一時(shí)刻的第1～第6步預(yù)測值，同時(shí)將原始數(shù)據(jù)中100～196 s的數(shù)據(jù)作為100～190 s中每一時(shí)刻的第1～第6步預(yù)測值的對照數(shù)據(jù)。

圖2 入侵機(jī)三維航跡圖Fig.2 3D track of the intruder

圖3 入侵機(jī)X軸位置變化Fig.3 X-axis position change of the intruder

圖4 入侵機(jī)Y軸位置變化Fig.4 Y-axis position change of the intruder

圖5 入侵機(jī)Z軸位置變化Fig.5 Z-axis position change of the intruder

4.2 參數(shù)設(shè)置

改進(jìn)滑動(dòng)窗多項(xiàng)式擬合法的階數(shù)m與滑動(dòng)窗長度n，均根據(jù)對航空器運(yùn)動(dòng)模式的判定自適應(yīng)取值，令T=1 s，q=6，依據(jù)航空器一般性能[15]和多項(xiàng)式預(yù)測模型使用經(jīng)驗(yàn)，對λ取值，設(shè)定λ=0.05 rad。

根據(jù)傳統(tǒng)滑動(dòng)窗多項(xiàng)式擬合法參數(shù)設(shè)置的不同，劃分四個(gè)對照組，為了對比兩個(gè)改進(jìn)措施對多項(xiàng)式擬合模型的影響，設(shè)定在傳統(tǒng)滑動(dòng)窗多項(xiàng)式擬合法中僅加入?yún)?shù)自適應(yīng)的改進(jìn)措施為對照組1。根據(jù)經(jīng)驗(yàn)，將基于傳統(tǒng)滑動(dòng)窗多項(xiàng)式擬合法的其余對照組參數(shù)設(shè)置為：對照組2，歷史數(shù)據(jù)個(gè)數(shù)n=2，多項(xiàng)式次數(shù)m=1；對照組3，n=3，m=1；對照組4，n=3，m=2。

4.3 誤差指標(biāo)

采用平均絕對誤差(mrerr)對航跡預(yù)測結(jié)果進(jìn)行分析，統(tǒng)計(jì)所有預(yù)測時(shí)間中每一步預(yù)測的絕對誤差并取平均值，能直觀地表現(xiàn)預(yù)測模型的預(yù)測性能，計(jì)算公式為

(10)

4.4 結(jié)果分析

首先，采用改進(jìn)的滑動(dòng)窗多項(xiàng)式擬合法，計(jì)算航跡數(shù)據(jù)，根據(jù)航跡預(yù)測實(shí)現(xiàn)步驟及參數(shù)設(shè)置，得到航空器在100～190 s內(nèi)每一時(shí)刻的第1～第6步預(yù)測值，并將預(yù)測值與原數(shù)據(jù)進(jìn)行比較，計(jì)算出絕對誤差，并繪制相應(yīng)的曲線；然后，采用傳統(tǒng)滑動(dòng)窗多項(xiàng)式擬合法，對航跡數(shù)據(jù)進(jìn)行計(jì)算，并設(shè)置參數(shù)，得到四組航空器在100～190 s內(nèi)每一時(shí)刻的第1～第6步預(yù)測值，計(jì)算每個(gè)對照組的預(yù)測值與原數(shù)據(jù)的絕對誤差，如圖6～圖10所示，并計(jì)算各組中每步預(yù)測值的平均絕對誤差，如表1所示。

圖6 改進(jìn)多項(xiàng)式擬合法航跡預(yù)測絕對誤差Fig.6 Absolute error of track predicted by improved polynomial fitting method

圖7 第一對照組航跡預(yù)測絕對誤差Fig.7 Absolute error of track prediction for first group

圖8 第二對照組航跡預(yù)測絕對誤差Fig.8 Absolute error of track prediction for second group

圖9 第三對照組航跡預(yù)測絕對誤差Fig.9 Absolute error of track prediction for third group

圖10 第四對照組航跡預(yù)測絕對誤差Fig.10 Absolute error of track prediction for fourth group

表1 平均絕對誤差對比

從圖6～圖10可以看出：平均絕對誤差值在減小，這是由于所加入的噪聲方差正比于兩機(jī)距離。綜合分析圖6～圖10與表1，可以得出：

(1) 本文改進(jìn)滑動(dòng)多項(xiàng)式擬合法應(yīng)用于航跡預(yù)測時(shí)，其預(yù)測精度要優(yōu)于傳統(tǒng)滑動(dòng)窗多項(xiàng)式擬合法。

(2) 對照組2在第1步和第2步的平均絕對誤差與本文改進(jìn)滑動(dòng)窗多項(xiàng)式擬合法的mrerr基本相同，但隨著預(yù)測步數(shù)的增加，對照組2中的兩個(gè)誤差指標(biāo)值具有比改進(jìn)滑動(dòng)窗多項(xiàng)式擬合法更大的增加量，表明本文改進(jìn)方法比傳統(tǒng)方法預(yù)測性能更加穩(wěn)定。

(3) 對照組1的mrerr要小于對照組2的，表明在傳統(tǒng)滑動(dòng)窗多項(xiàng)式擬合法中自適應(yīng)調(diào)整多項(xiàng)式參數(shù)會(huì)得到精度更高的預(yù)測值，圖7～圖8也證明了此結(jié)論的正確性；通過對照組1與改進(jìn)多項(xiàng)式擬合法的mrerr指標(biāo)值對比，證明了為每個(gè)預(yù)測值提供更合適的多項(xiàng)式擬合方程的改進(jìn)措施，提高了航跡預(yù)測的精度，圖6～圖7也證明了此結(jié)論的正確性。

(4) 對照組3中的每步mrerr值相對于前一步mrerr值的增長序列與對照組1大體相同，這是由于對照組3與對照組2多項(xiàng)式階數(shù)相同所致，但由于滑動(dòng)窗長度更大，對照組2中的預(yù)測誤差均值遠(yuǎn)大于對照組1。這也印證了式(7)的正確性。

5 結(jié) 論

(1) 在使用多項(xiàng)式擬合模型同時(shí)預(yù)測多個(gè)連續(xù)未來值時(shí)，依據(jù)預(yù)測值與當(dāng)前值的時(shí)間間隔為預(yù)測值挑選歷史數(shù)據(jù)組，為各預(yù)測值構(gòu)造合適的多項(xiàng)式擬合方程，可提高航跡預(yù)測精度。

(2) 在單個(gè)維度中目標(biāo)航跡點(diǎn)所構(gòu)成的時(shí)間-位置折線圖上，依據(jù)相鄰折線所形成角度大小及角度的連續(xù)變化趨勢進(jìn)行運(yùn)動(dòng)模式的判斷，實(shí)現(xiàn)擬合多項(xiàng)式階數(shù)自適應(yīng)與滑動(dòng)窗長度自適應(yīng)，這樣預(yù)測的三維航跡精度更高。

(3) 改進(jìn)滑動(dòng)多項(xiàng)式擬合法的預(yù)測精度優(yōu)于傳統(tǒng)滑動(dòng)窗多項(xiàng)式擬合法，且進(jìn)行多步預(yù)測時(shí)預(yù)測性能也更加穩(wěn)定。未來的研究工作將更多聚焦于本文方法所預(yù)測4D航跡在無人機(jī)飛行防相撞實(shí)踐中的驗(yàn)證。