鄭天宇,姚 郁,賀風華
(哈爾濱工業大學 航天學院,哈爾濱 150080)
臨近空間高超聲速目標因具有飛行空域廣、機動能力強和全球到達的特性,已成為國防安全的新型威脅. 此類目標多采用乘波體、升力體、翼面融合體等氣動構型,具有極強的機動能力,航跡設計具有較大的自由度,可進行多種類型的規避和突防策略穿越已知防區,甚至可通過大范圍的側向機動更換擬攻擊目標,給航跡估計帶來了新的挑戰.
近年來,機動目標的航跡估計問題得到了廣泛而深入的研究.研究者們根據不同的目標運動特性,設計相應的機動模型,并提出了各類估計方法[1-4].文獻[1]給出了估計中常用的模型,包括恒速(constant velocity,CV)模型、恒加速度(constant acceleration,CA)模型、恒轉彎速率(constant turn rate,CTR)模型、非對稱法相加速度模型、Singer模型、當前統計模型等.當目標的動力學模型為線性且量測噪聲為零均值高斯噪聲時,卡爾曼濾波(Kalman filter,KF)可獲得最優的航跡估計.但是在實際估計問題中,目標的動力學和機動模型往往是非線性的,因此研究者們提出了一系列基于非線性濾波的估計方法[3].針對機動目標航跡估計問題,文獻[5]提出了一種基于改進Sage-Husa自適應卡爾曼濾波的機動目標跟蹤方法,通過設計判斷準則調整自適應濾波參數,抑制了目標跟蹤誤差的積累,提高了穩態濾波精度.文獻[6]結合Singer機動模型與模糊推理,提出了一種多參數融合的Singer模型并將其用于目標機動估計,取得了較高的估計精度.文獻[7]提出了一種帶有加速度預測的非線性估計算法,通過最小二乘估計器對滑動窗口內目標機動的檢測,以一定的策略調整目標機動方差的先驗信息,提升了機動估計的精度.文獻[8]提出了一種變速率粒子濾波算法,可遞歸的更新狀態序列的后驗分布,從而解決了變維狀態推斷問題,可用于機動模型不確定的目標跟蹤問題.文獻[9]將成本回溯自適應輸入和狀態估計(retrospective-cost-based adaptive input and state estimation,RCAISE)用于機動目標的跟蹤,可通過優化回溯性能使得估計的機動加速度逼近目標真實的加速度.針對強機動目標跟蹤時的濾波增益滯后和過程噪聲不精確的問題,文獻[10]提出了一種改進的自適應無跡卡爾曼濾波方法,通過在每步估計中加入自適應的比例因子更新過程噪聲協方差陣提升了估計精度和動態性能.針對機動目標跟蹤時的模型失配問題,文獻[11]提出了一種自適應高階無跡卡爾曼濾波方法,通過引入一種由預測殘差估計協方差矩陣求取最優自適應因子減小了動態模型誤差對估計的影響.這些方法均基于卡爾曼濾波或粒子濾波框架提出,使用單一的模型對目標航跡進行估計,并通過UT變換、自適應等技術應對機動模型的不確定性.對于臨近空間高超聲速目標的航跡估計問題,單一機動模型難以準確描述目標復雜的機動.在應用這些方法時,目標機動將帶來大范圍的模型不確定性,常見自適應機制難以應對,進而導致估計精度下降.但是,目標航跡是“設計-優化-控制”的結果,統計意義上存在一定的規律,可通過引入循環神經網絡(recurrent neural networks,RNNs)對目標機動機動特性進行識別.
循環神經網絡在20世紀80年代由Lipton等[12]提出.近年來在語音識別、文本生成、機器翻譯等序列數據的處理領域得到了廣泛的應用[13].目前,循環神經網絡與非線性濾波的聯合使用方面已有一些研究成果,提出了一系列循環神經網絡與擴展卡爾曼濾波(extended Kalman filter, EKF)、無跡卡爾曼濾波(unscented Kalman filter,UKF)、粒子濾波(particle filter,PF)等結合的估計方法[14-16].文獻[14]考慮視頻有損壓縮-解碼中的偽影抑制問題,設計了一種深度卡爾曼濾波器,該方法使用兩個卷積神經網絡替代了卡爾曼濾波的狀態和誤差協方差陣預測步驟,提升了視頻解碼的質量.文獻[15]提出了一種基于KF和循環神經網絡的人體姿勢估計算法,通過使用KF融合兩個循環神經網絡的跟蹤結果,提升了跟蹤的位置和速度精度.文獻[16]提出了一種用于人體姿勢估計的KF-LSTM網絡結構,該方法將3個長短時記憶(long short-term memory,LSTM)網絡嵌入KF中并通過訓練學習運動和噪聲模型,實現了優于現有結果的姿勢估計精度.考慮機器人的視覺定位問題,文獻[17]通過將系統模型和粒子濾波算法嵌入神經網絡提出了PF-net網絡結構,該方法由序列數據端對端的學習和優化粒子濾波中使用的模型.文獻[18]考慮噪聲條件下的手勢識別問題,通過將UKF嵌入LSTM網絡設計了NIUKF-LSTM網絡,可有效抑制量測噪聲提升識別精度.這些研究中所涉及的動力學模型均較為簡單,濾波方法中的模型部分往往被循環神經網絡完全的替代并通過訓練學習.但是,臨近空間高超聲速目標的動力學模型存在較為復雜的非線性特性,難以直接通過訓練學習.因此,需尋求新的RNN-EKF結合方式.
本文考慮臨近空間高超聲速目標的航跡估計問題,提出基于可學習擴展卡爾曼濾波(Learnable EKF,L-EKF)的航跡估計方法.通過將循環神經網絡與擴展卡爾曼濾波方法深度嵌合,識別目標復雜的機動特征,提升濾波方法對復雜目標機動的應對能力,實現高精度的航跡估計.
臨近空間高超聲速目標可在臨近空間高度長時間滑翔飛行,具有非慣性的航跡形式和復雜的策略性機動,常見的CA、CV、Singer模型難以準確描述其機動特征.因此,本文對其機動特性進行分析和描述:首先,給出臨近空間高超聲速目標的動力學模型;然后,分析其機動特性,建立參數化的航跡估計模型,進而給出航跡估計問題的數學描述.
忽略地球自轉,高超聲速目標的動力學模型可表示為[19]:

(1)
(2)
(3)
(4)
(5)
(6)
式中:r為目標到地心的距離;θ為經度;φ為緯度;v為速度;γ為航跡傾角;ψ為航跡偏角;σ為傾側角;m為目標質量;Sref為參考面積;g為重力加速度;ρ=ρ0e-βh為大氣密度;h=r-Re為目標高度,Re=6 378 135 m為地球半徑;ρ0為海平面處大氣密度;β為大氣密度系數;CL、CD分別為升力、阻力系數.
由于本文所考慮的目標是非合作的,其精確的升力和阻力系數在實際的航跡估計中難以獲得,因此使用馬赫數無關模型[20]描述其升力和阻力,表示如下:
式中CD0、K分別為馬赫數無關的氣動參數.
目標的升阻比定義為
(7)

歸一化升力系數定義為
則目標的升力系數CL和阻力系數CD可表示為:
(8)
(9)
聯立式(1)~(6), 式(8)、(9)可得高超聲速目標的動力學模型:
(10)
在航跡估計過程中,歸一化升力系數cl和傾側角σ往往難以直接獲得.同時,由于目標具有復雜且不確定的機動形式,cl和σ存在復雜時變特性,難以被直接識別.換而言之,動力學模型(10)無法在航跡估計中直接使用.因此,本文針對目標的機動特性機型分析,選擇可識別的參數替代cl和σ,建立參數化的目標機動模型.


(11)
將式(8)、式(11)代入式(5)可得:
(12)

(13)
由式(13)可知, 準平衡滑翔目標h-v之間的關系由clcosσ決定.因此,選擇clcosσ作為目標航跡的第1個特征參數,記作:
λ1=clcosσ.
(14)
對式(12)求關于時間的導數,可得
(15)
將式(15)代入式(1),可得
(16)
(17)
由于參數λ1,λ2僅與cos(σ)有關,無法表征目標側向機動的方向,本文將傾側角σ的符號選為目標航跡的第3個特征參數,即有
(18)
目標飛行過程中,其機動方向往往變化緩慢.為保證模型連續并描述其變化特性,將λ3近似為tanh函數,并使用一階馬爾科夫模型表述,即有:
λ3≈tanh(kχχ),
(19)
(20)
式中:kχ為側向機動方向模型的增益系數,可將其設為一個很大的正數;Tχ為側向機動方向模型的時間常數.
綜合式(10)、(14)、(17)~(20)可得參數化的目標機動模型:

(21)
由式(21)可知,通過參數λ1,λ2和參數λ3的時間常數Tχ即可描述目標的機動特性.這3個參數與目標的高度、速度、機動周期等特性具有較強的關聯性,易于使用循環神經網絡根據目標量測信息直接識別.
下面,給出航跡估計問題的標準形式.由于目標只有位置信息可被直接測量,因此選擇狀態為x=[r,θ,φ,v,γ,ψ,χ]T,量測為目標位置z=[r,θ,φ]T,輸入為需識別的3個特征參數u=[λ1,λ2,Tχ]T,目標的航跡估計的標準形式為

(22)

實現臨近空間高超聲速目標航跡估計的關鍵在于處理非慣性的航跡形式和復雜的策略性機動. 臨近空間高超聲速目標機動特性分析與描述中,本文已經建立了可描述目標航跡特性的參數化機動模型(22),在航跡估計中需對這些特征參數進行識別.此外,受臨近空間復雜氣動環境的影響,目標很多的高階動力學特性無法被模型描述,在航跡估計中也需對這些未建模的動力學特性進行處理.因此,本文考慮臨近空間高超聲速目標的航跡估計問題,提出可學習擴展卡爾曼濾波(Learnable extended Kalman filter,L-EKF)方法.首先,分析擴展卡爾曼濾波(EKF)方法在應對大范圍模型和參數不確定性時存在的問題;然后,分別設計輸入修飾網絡(input modification network,IMN)和增益修飾網絡(gain modification network,GMN),并將其嵌入EKF中,用以識別和補償估計中存在的各類不確定性;最后,給出完整的可學習擴展卡爾曼濾波算法.
考慮如下離散非線性系統:
式中:xk為狀態;zk為量測;uk為輸入;wk、vk分別為過程噪聲和量測噪聲,均為零均值的高斯噪聲.
則經典的擴展卡爾曼濾波算法[21]如下(如圖1所示).

圖1 擴展卡爾曼濾波
Step1預測狀態和誤差協方差陣:
(23)
(24)

Step2計算次優卡爾曼增益:
(25)
(26)
(27)

Step3更新狀態和誤差協方差矩陣:
Pk|k=(I-KkHk)Pk|k-1.
在實際的航跡估計中,通常情況下EKF不是最優的估計器,其原因主要來源于模型的不確定性和線性化誤差,具體體現在以下3個方面:
1)模型的未建模不確定性. 在實際的航跡估計中,目標很多的高價動力學特性無法被精確建模,即估計模型存在未建模不確定性.此類不確定性帶來的過程噪聲通常是非高斯的,難以被過程噪聲協方差陣Qk精確的描述.
2)模型的輸入不確定性.如式(22)所述,在進行航跡估計時模型需要參數輸入uk.由于所考慮的目標是非合作的,實際航跡估計中,這些參數難以被精確獲得,即uk存在不確定性.
3)模型的線性化誤差. 在EKF中,每一步估計均需對模型進行工作點處線性化,由此引入了模型線性化誤差.
這些不確定性和誤差在EKF框架下無法被準確的描述和處理,將影響航跡估計的性能.但從更為廣義的角度,它們也是目標特性的一部分,對于某類具體目標而言存在一定的統計特性,可使用已有數據訓練循環神經網絡進行識別和補償.因此,本文設計兩個循環神經網絡并將其嵌入EKF中:一個為輸入修飾網絡,用于應對模型的輸入不確定性;另一個為增益修飾網絡,用于應對模型的未建模不確定性和線性化誤差.

式中κk為IMN的隱狀態.
為了應對目標的周期性機動,IMN需要長時記憶能力.因此,本文將IMN設計為帶有批規范化的LSTM編碼網絡結構. 如圖2所示,網絡由2個全連接層(Dense),2個LSTM[22]層,1個輸出層(Output)和1個加權層(Weight)串聯而成,并在每個全連接層后加入批規范化(batch normalization,BN)[23]層.全連接層用于在維度上豐富輸入的特征信息,記作:
δk,i=Densei(μk,i)=tanh(ωiμk,i+bi).
式中:Densei為第i個全連接層;μk,i為Densei層的輸入;δk,i為Densei層的隱狀態,也是輸出;ωi、bi分別為Densei層的可訓練參數.

圖2 輸入修飾網絡
LSTM層用于在時間線上對特征信息進行“編碼”,以表征各時刻輸入時間上的相關性,記作:
[δk,i,LSTM,Ck,i,LSTM]=
LSTMi(δk-1,i,LSTM,Ck-1,i,LSTM,μk,i,LSTM).
式中:LSTMi為第i個LSTM層;μk,i,LSTM為LSTMi的輸入;δk,i,LSTM為LSTMi的隱狀態,也是輸出;Ck,i,LSTM為LSTMi的胞狀態(cell states).
批規范化層則用于抑制過擬合,提升網絡的泛化能力,其數學表述為

輸出層用于解算uk的加權系數,需根據模型輸入uk的特性設計其值域.綜上所述由模型(22)可知u=[λ1,λ2,Tχ]T,這3個參數均為正數且變化幅度不大,因此輸出層使用帶有偏置的tanh激活函數,表示為
κi,k=Outputimn(μk)=
加權層用于根據之前解算的加權系數計算修飾后的模型輸入,使用乘性加權方式,設計為
式中*為按元素乘法.
綜上所述,輸入修飾網絡的數學表達如下:

式中Gk為GMN的隱狀態.
類似于IMN,本文將GMN設計為帶有批規范化的LSTM編碼網絡結構. 如圖3所示,網絡由2個全連接層(Dense),2個LSTM層,1個輸出層(Output)和1個加權層(Weight)串聯而成,并在每個全連接層后加入批規范化層(BN). 全連接層用于在維度上豐富輸入的特征信息;LSTM層用于在時間線上對特征信息進行“編碼”,以表征每一時刻特征間的時間相關性;批規范化層用于應對過擬合,提升網絡的泛化能力.
輸出層用于解算次優卡爾曼增益的加權系數,由于模型(22)中各狀態量級差距較大,為保證加權系數具有足夠的權重,使用指數函數作為輸出層的激活函數,則有
κg,k=Outputgmn(μk)=


圖3 增益修飾網絡
加權層用于根據之前解算的加權系數計算修飾后的次優卡爾曼增益,使用乘性加權方式,設計為

綜上所述,增益修飾網絡的數學表達如下:
下面給出可學習擴展卡爾曼濾波(L-EKF)算法,如圖4所示,本文在EKF的狀態和誤差協方差預測之前嵌入了IMN,在EKF的狀態和誤差協方差更新之前嵌入了GMN.由于L-EKF對過程和量測噪聲協方差陣Qk和Rk不敏感,本文在算法設計中用定常的噪聲協方差陣Q,R代替了Qk,Rk,以簡化參數設計.完整的L-EKF算法如下.
Step1輸入修飾:
Step2預測狀態和誤差協方差陣:

Step3計算次優卡爾曼增益:

Step4增益修飾:
Step5更新狀態和誤差協方差矩陣:
注釋1L-EKF方法給出了一種新的RNN-EKF聯合使用方式.在各類不同的估計問題中,IMN和GMN可設計為任意結構的循環神經網絡,包括但不限于Simple RNN, LSTM, GRU以及基于這些結構的自動編碼網絡.具體的網絡結構需根據估計問題的特點和復雜程度設計.

圖4 可學習擴展卡爾曼濾波
本文將所提出的L-EKF方法應用于臨近空間高超聲速目標的航跡估計.首先,對L-EKF方法進行訓練;然后,通過幾個仿真案例將L-EKF方法與擴展卡爾曼濾波(EKF)、自適應擴展卡爾曼濾波(Adaptive EKF,AEKF)等航跡估計方法進行對比,分析其航跡估計性能.
L-EKF是一種“端對端”方法,訓練中不需作為對數據進行標注,直接使用真實航跡的狀態x和輸入u作為標稱值.本文將一臺裝有Intel Core i7-8700K CPU,Nvidia TITAN Xp顯卡,32 GB內存的計算機,并基于Python 3.7+Tensorflow+Keras平臺對網絡進行訓練.訓練的具體設置如下.
3.1.1 網絡參數
LSTM2層的神經元數分別為256-32-256-256;GMN的結構如圖3所示,Dense1,Dense2,LSTM1,LSTM2層的神經元數分別為128-16-128-128.
3.1.2 數據集
使用一個由1 080條航跡構成的數據集進行訓練.數據集使用CAV-H[24]的氣動數據生成,包含不同cl的準平衡滑翔航跡,并考慮CV、CA、轉彎、方波等多種類型的側向機動.隨機選擇80%的航跡作為訓練集,10%的航跡作為驗證集,10%的航跡作為測試集.雖然可使用任意長度的航跡數據序列訓練L-EKF,但受算力所限將訓練中使用的航跡切割為長度100的定常片段.
3.1.3 損失函數
損失函數為平衡各狀態的量級,使用帶加權的平均絕對誤差損失函數(weighted mean absolute error,WMAE),設計如下:
式中:wx=[50,1 000,1 000,0.5,10 000,5 000],wu=[1,1,0.01],c=0.5為加權系數;xi、ui分別為航跡狀態和輸入的標稱值.
3.1.4 優化器
使用RMSprop[25]優化器,訓練20 epochs,學習率(learning rate)設置為0.001.
3.1.5 訓練結果
訓練結果如圖5所示,可見經過充分的訓練L-EKF的WMAE損失函數收斂.訓練集和測試集上的損失函數基本相等,說明所提出的方法具有良好的泛化能力,可應對未知的目標機動.

圖5 訓練結果
本文通過3個典型的目標機動場景,將所提出的L-EKF方法與傳統的EKF、AEKF方法對比,分析其航跡估計性能.假設僅有目標的高度h、經度θ、緯度φ信息可通過預警系統量測獲得,且量測信息帶有均值為0,均方根誤差分別為2 000 m, 0.01 rad, 0.01 rad 的高斯白噪聲.航跡仿真條件見表1,場景設定如下.

表1 仿真條件
場景1縱向做變歸一化升力系數的準平衡滑翔;側向701~1 300 s做過載為3 g的CA機動,其余時間無機動.
場景2縱向以恒定歸一化升力系數做準平衡滑翔;側向0~1 200 s做過載為1 g的方波機動,1 201~1 800 s做過載為5 g的大范圍規避機動,其余時間無機動.
場景3縱向做變歸一化升力系數的準平衡滑翔;側向401~800 s做過載為3 g的CA機動,1 501 s之后做過載為5 g的大范圍規避機動,其余時間無機動.
仿真結果見表2和圖6~8所示.由表2可知,在所有仿真場景中,L-EKF方法具有比EKF和AEKF更高的估計精度.由圖6~8可知,當目標不機動時,3種方法均有較高的估計精度.當目標機動時,受模型不確定性的影響,EKF的估計精度有著較大的下降,L-EKF和AEKF方法均可以一定程度上適應這種不確定性.由于IMN和GMN可準確的識別目標的機動特性并補償由其引發的估計誤差,L-EKF方法的估計精度和動態性能均優于AEKF方法.運算速度方面,由于嵌入了兩個循環神經網絡,L-EKF方法略慢于EKF和AEKF方法,但每步估計的計算耗時仍處于同一量級,說明L-EKF方法具有較好的實時性.綜上所述,所提出的L-EKF方法可準確識別未知的模型參數并補償由其帶來的估計誤差,有效應對臨近空間高超聲速目標復雜的機動,具有更高的估計精度和更好的估計動態性能.

表2 估計結果

圖6 場景1仿真結果

圖7 場景2仿真結果

圖8 場景3仿真結果
1)分析了臨近空間高超聲速目標的機動特性,建立了參數化的目標機動模型,解決了經典機動模型難以描述高超聲速目標運動特性的問題;
2)針對目標復雜機動導致的參數和模型不確定性,提出了可學習擴展卡爾曼濾波(L-EKF)方法.通過設計和訓練兩個循環神經網絡(IMN和GMN),并將其嵌入EKF中,實現了對目標復雜機動引起的參數和模型不確定性的在線識別與動態補償;
3)通過典型機動條件下的對比分析,可知所提出的L-EKF方法可有效應對目標復雜的機動,具有比傳統方法更高的估計精度和更優的估計動態性能.