張建祥,甘旭升,周志靖, 楊捷
(1.西京學院 理學院,西安 710123) (2.空軍工程大學 空管領航學院,西安 710051)
航跡是指航空器在飛行過程中按時間順序所經歷的全部或部分空間位置點的集合,將時間與各個空間點位一一對應后的組合被稱為4D航跡。航跡預測是指根據航空器初始飛行計劃、歷史航跡信息與經驗信息,采用相應預測方法,對航空器未來時間段內將要產生的航跡點進行預測。高精度的4D航跡預測是無人機飛行防相撞空中威脅態勢預警技術的基礎,是提高空中交通管理效率與安全水平的重要技術手段。
卡爾曼濾波是由 R.E.Kalman[1]于1960年推導出來的一種最小方差最優估計算法,主要用于導航系統多傳感器信息融合處理和目標狀態估計。1998年,L.Tong等[2]首次應用卡爾曼濾波算法解決航跡預測問題,對相鄰時刻目標位置及其瞬時速度進行預測,從一定程度上改善了預測精度,但在實際應用中存在非線性系統的適用性問題,預測性能還有較大提升空間。之后,又有諸多研究人員先后提出改進型的基于殘差均值的交互式多模型跟蹤濾波算法處理航跡預測問題,該濾波算法對機動目標追蹤預測有較高精度,但需要做大量的數學建模工作,例如航空器運動模型、運行環境模型等,并最終融合各種模型,形成航跡預測模型,工作量非常巨大[3-6]。此外,引入數據挖掘方法對航空器歷史航跡數據進行處理、分析,挖掘數據中隱含的模式,從而預測航空器未來位置,并已在航空器4D航跡預測領域得到應用,比較成熟的方法,例如神經網絡算法,存在局部收斂問題和網絡結構確定問題。美國艾姆斯研究中心與美國管制員協會聯合研發了中央終端區管制自動化系統[7],其航跡預估模塊是核心部分,可生成空域內所有航空器未來30 min內的航跡。PATs[8]是歐洲開發的航跡預測工具,可預先評估所發出的管制指令對航空器的影響,提高歐洲地區空中交通管理效率與安全水平。國內的航空器4D航跡預測研究起步較晚。張軍峰等[9]提出了基于連續動態模型與離散動態模型的4D航跡預測方法,并將其用于離場航跡預測,但該方法復雜度較高,無法滿足實時性需求;王建忠等[10]根據點融合進近的特點,提出航空器4D航跡預測方法,并采用遺傳算法進行求解,但該方法僅能處理簡單問題,無法預測多跑道航跡的復雜情況;張振興等[11]提出基于實時反饋長短期記憶神經網絡的4D航跡預測模型,但該方法確定網絡參數存在較大難度,所建模型通常無法保證精度。
本文提出采用滑動窗多項式擬合法解決航空器4D航跡預測問題。對滑動窗多項式擬合法進行改進,構建4D航跡預測模型,力求給出一種精度高、簡單易行且計算簡便的航跡預測方法,并通過實例驗證其有效性和可行性。

曲線擬合的具體作法:針對給定數據點{Xi|i=0,1,…,n},在給定的函數類Φ中,求f(ti)∈Φ,使誤差ri=f(ti)-Xi(i=0,1,…,n)的平方和最小,即
(1)
從幾何意義上講,就是找尋出一條曲線y=f(t),使之與給定點{Xi|i=0,1,…,n}的距離平方和為最小。函數f(t)稱為擬合函數或最小二乘解,求擬合函數f(t)的方法稱為曲線擬合的最小二乘法。

(2)

A=(PTP)-1PTX
(3)

易證式(3)中(PTP)為對稱正定矩陣,故存在唯一解。求解式(3)可得擬合多項式
(4)
由式(4)可對(tn+dt)時刻的目標航跡進行預測,即采用多項式擬合法預測航跡:
(5)
固定式(5)中的n值,即預測模型取固定數量的數據點,并隨著周期更新,不斷接收數據點,剔除相同數量的舊數據點,即在多項式擬合法中應用滑動窗思想,模型將具有動態預測能力。
為了最大程度上利用歷史數據信息,提高算法預測性能,對傳統滑動窗多項式擬合法提出如下改進措施。
傳統滑動窗多項式擬合法的優點在于采用遞推法最大程度地利用有限的歷史數據信息,而本文對滑動窗多項式擬合法的改進之處則是進一步強化了該優點。傳統滑動窗多項式擬合法基于歷史數據組構建了一個多項式擬合方程,依據該多項式擬合方程計算所有預測值,但每個預測值與歷史數據組的時間距離均不同,僅僅依靠同一個多項式擬合方程,并不能完全表達出每個未來值受歷史數據組的影響程度。本文在同時預測多個連續未來值時,為每個預測值構造了更為合適的多項式擬合方程,以便更準確地體現歷史數據組對每個未來值的不同影響,將得到比傳統滑動窗多項式擬合法更高的預測精度。
為每一個預測值構建更合適的多項式擬合方程的關鍵在于挑選合適的歷史數據組,本文基于每個預測值與當前值的時間距離來挑選歷史數據組,充分考慮多項式擬合法的時間序列性質。在實際應用中,假設當前值為Xk,滑動窗寬度為n,即預測模型中的歷史數據個數為n,預測當前值之后1,2,…,q時刻的未來值,針對當前值之后第s(1≤s≤q)時刻的未來值進行預測時,其與當前值的時間距離為s,則其歷史數據組的選擇應為{Xk-s(n-1),Xk-s(n-2),…,Xk-s,Xk},再依據式(5),即可為每個預測值提供更合適的多項式擬合方程。
本文所提出的滑動窗多項式參數自適應主要是指擬合多項式階數與滑動窗長度的自適應。
2.2.1 擬合多項式的階數自適應
根據多項式的最小二乘擬合原理可知,在計算預測模型時,首先要確定擬合多項式的階數m。要獲得高精度的航跡預測值,必須判明目標的運動狀態,使擬合多項式的階數與其相適應。本文將目標運動模式劃分為類直線運動與曲線運動,設定目標處于類直線運動模式與曲線運動模式時,對應的擬合多項式階數m分別為1和2。
對于三維空間中目標的運動模式是否存在類直線與曲線之間的突變以及如何進行突變的判定問題,將其簡化到三個單維空間之中進行解決,根據目標在單個維度之中的航跡點繪制時間-位置折線圖,取相鄰折線所形成角度的大小以及角度的連續變化趨勢為判定依據。
在目標航跡序列中取X軸上相鄰的四個航跡點為xk-3,xk-2,xk-1,xk,由于具有相同的時間間隔,根據上述四個航跡點可直接繪制成時間-位置折線圖,如圖1所示,α1為目標在X軸上航跡點xk-3,xk-2,xk-1形成相鄰折線的夾角,α2為目標在X軸上航跡點xk-2,xk-1,xk形成相鄰折線的夾角。
根據向量夾角公式可得出圖1中兩個連續夾角的計算方程為[13]
(6)
式中:arccos為三角函數中反余弦符號;‖‖2為向量的2-范數符號[14]。

圖1 時間-位置折線圖Fig.1 Time-position line chart
可通過設定一個參數λ來進行目標運動模式與相關突變情況的判定:
若目標航跡滿足α1>λ,α2>λ,則目標在X軸上處于曲線運動模式,擬合多項式階數取m=2;若目標航跡滿足α1>λ,α2≤λ,則目標在X軸上是從曲線運動模式轉為類直線運動模式,擬合多項式階數取m=1;若目標航跡滿足α1≤λ,α2>λ,則目標在X軸上是從類直線運動模式轉為曲線運動模式,擬合多項式階數取m=2;若目標航跡滿足α1<λ,α2<λ,則目標在X軸上是處于類直線運動模式,擬合多項式階數取m=1。
目標在Y軸與Z軸上的運動模式與相關突變情況的判定,與目標在X軸上的判定原理相同,此處不再贅述。
2.2.2 多項式滑動窗長度自適應
因為預測是基于有限的歷史數據進行的,所以必須選擇所使用時間窗寬度。如果時間窗寬度較大,即歷史數據太多,則預測算法的運算量加大,導致預測緩慢,且如果目標運動軌跡變化較快,將導致實際目標運動軌跡與時間多項式函數之間的擬合效果惡化,則預測偏差增加;如果時間窗寬度太小,則數據量過少,所計算出的多項式函數將不能描述目標的運動特性,預測效果也不理想。因此,確定滑動窗寬度時,一方面要考慮不能使預測模型計算太多的歷史數據,即滑動窗寬度不可太寬,另一方面則要兼顧目標運動軌跡變化的快慢,即須考慮所選取的擬合多項式階數。
綜上所述,為了進一步減少擬合誤差,提高預測精度,依據多項式擬合原理得出擬合多項式滑動窗長度n與階數m之間的關系式(式(7)),而且在低階數多項式擬合法的位置預測研究中,李亞寧[12]已通過仿真證明了其合理性。
n=m+1
(7)
采用改進的滑動窗多項式擬合法對航空器航跡進行預測,其過程是先在直角坐標系中X,Y,Z軸上分別對航空器航跡進行預測,而后形成對航空器三維空間航跡的總體預測。具體步驟如下:
(1) 按照改進方法在X,Y,Z軸上為每個預測值挑選合適的歷史數據組。
(2) 依據式(6)及目標運動模式與相關突變情況的判定結果,確定X,Y,Z軸上的擬合多項式階數mx,my,mz,滑動窗長度nx,ny,nz。航空器在X,Y,Z軸上的運動軌跡多項式擬合函數為
(8)

(9)
式中:T為更新周期。
求得式(9)中X,Y,Z軸上航空器航跡預測值,即可得航空器三維空間中航跡總體預測值。
(4) 采用遞推算法,重新為每個預測值挑選合適的歷史數據組,并利用滑動窗口不斷更新最新的X,Y,Z軸上航空器航跡數據,重新確定X,Y,Z軸上的擬合多項式階數mx,my,mz與滑動窗寬度nx,ny,nz,求得基于線性最小方差估計的多項式預測模型,預測下一周期航空器的航跡。
在MATLAB 2014a計算環境中,分析改進滑動窗多項式擬合法預測航跡的性能,并與傳統滑動窗多項式擬合法進行對比。由于后者無法自適應調整多項式階數與滑動窗長度,在仿真分析中,將根據階數與滑動窗長度不同,設置多個對照組。
實驗使用數據來自于通過六階貝塞爾曲線法構建的入侵機航跡數據,即入侵機相對無人機的航跡數據,如圖2~圖5所示。在所有的原始數據中,將其中0~99 s的數據只作為歷史數據使用,預測航跡為100~196 s中每一時刻的第1~第6步預測值,同時將原始數據中100~196 s的數據作為100~190 s中每一時刻的第1~第6步預測值的對照數據。

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

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

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

圖5 入侵機Z軸位置變化Fig.5 Z-axis position change of the intruder
改進滑動窗多項式擬合法的階數m與滑動窗長度n,均根據對航空器運動模式的判定自適應取值,令T=1 s,q=6,依據航空器一般性能[15]和多項式預測模型使用經驗,對λ取值,設定λ=0.05 rad。
根據傳統滑動窗多項式擬合法參數設置的不同,劃分四個對照組,為了對比兩個改進措施對多項式擬合模型的影響,設定在傳統滑動窗多項式擬合法中僅加入參數自適應的改進措施為對照組1。根據經驗,將基于傳統滑動窗多項式擬合法的其余對照組參數設置為:對照組2,歷史數據個數n=2,多項式次數m=1;對照組3,n=3,m=1;對照組4,n=3,m=2。
采用平均絕對誤差(mrerr)對航跡預測結果進行分析,統計所有預測時間中每一步預測的絕對誤差并取平均值,能直觀地表現預測模型的預測性能,計算公式為
(10)

首先,采用改進的滑動窗多項式擬合法,計算航跡數據,根據航跡預測實現步驟及參數設置,得到航空器在100~190 s內每一時刻的第1~第6步預測值,并將預測值與原數據進行比較,計算出絕對誤差,并繪制相應的曲線;然后,采用傳統滑動窗多項式擬合法,對航跡數據進行計算,并設置參數,得到四組航空器在100~190 s內每一時刻的第1~第6步預測值,計算每個對照組的預測值與原數據的絕對誤差,如圖6~圖10所示,并計算各組中每步預測值的平均絕對誤差,如表1所示。

圖6 改進多項式擬合法航跡預測絕對誤差Fig.6 Absolute error of track predicted by improved polynomial fitting method

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

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

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

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

表1 平均絕對誤差對比
從圖6~圖10可以看出:平均絕對誤差值在減小,這是由于所加入的噪聲方差正比于兩機距離。綜合分析圖6~圖10與表1,可以得出:
(1) 本文改進滑動多項式擬合法應用于航跡預測時,其預測精度要優于傳統滑動窗多項式擬合法。
(2) 對照組2在第1步和第2步的平均絕對誤差與本文改進滑動窗多項式擬合法的mrerr基本相同,但隨著預測步數的增加,對照組2中的兩個誤差指標值具有比改進滑動窗多項式擬合法更大的增加量,表明本文改進方法比傳統方法預測性能更加穩定。
(3) 對照組1的mrerr要小于對照組2的,表明在傳統滑動窗多項式擬合法中自適應調整多項式參數會得到精度更高的預測值,圖7~圖8也證明了此結論的正確性;通過對照組1與改進多項式擬合法的mrerr指標值對比,證明了為每個預測值提供更合適的多項式擬合方程的改進措施,提高了航跡預測的精度,圖6~圖7也證明了此結論的正確性。
(4) 對照組3中的每步mrerr值相對于前一步mrerr值的增長序列與對照組1大體相同,這是由于對照組3與對照組2多項式階數相同所致,但由于滑動窗長度更大,對照組2中的預測誤差均值遠大于對照組1。這也印證了式(7)的正確性。
(1) 在使用多項式擬合模型同時預測多個連續未來值時,依據預測值與當前值的時間間隔為預測值挑選歷史數據組,為各預測值構造合適的多項式擬合方程,可提高航跡預測精度。
(2) 在單個維度中目標航跡點所構成的時間-位置折線圖上,依據相鄰折線所形成角度大小及角度的連續變化趨勢進行運動模式的判斷,實現擬合多項式階數自適應與滑動窗長度自適應,這樣預測的三維航跡精度更高。
(3) 改進滑動多項式擬合法的預測精度優于傳統滑動窗多項式擬合法,且進行多步預測時預測性能也更加穩定。未來的研究工作將更多聚焦于本文方法所預測4D航跡在無人機飛行防相撞實踐中的驗證。