張凱， 熊家軍， 付婷婷， 席秋實

(1. 空軍預警學院 研究生大隊， 武漢 430019； 2. 空軍預警學院 預警情報系， 武漢 430019)

作為新一代跨大氣層空天飛行器，高超聲速滑翔目標(Hypersonic Gliding Vehicle, HGV)結(jié)合了航天器與航空器的特征，具有高速高機動、高精度、大航程等特點，可實施遠距離機動快速打擊[1-2]。隨著美國分別于2011年和2017年成功試射陸基和潛射型高超聲速滑翔彈頭(AHW和IRCPS)，標志著HGV武器已初步具備作戰(zhàn)能力，這給中國空天安全帶來了巨大威脅。面對這類目標早期預警和攔截中的目標跟蹤問題，有必要結(jié)合HGV的動力學特性進一步研究相關(guān)理論。

不同于傳統(tǒng)高速目標，HGV沿非慣性彈道機動飛行，地基傳感器存在探測時間短、測量精度低等問題，使用傳統(tǒng)跟蹤算法難以對HGV實現(xiàn)準確跟蹤。從濾波器應用角度來看HGV跟蹤是典型的非線性濾波問題，其估計性能很大程度上取決于目標狀態(tài)模型的精細程度。一般來說，合理構(gòu)建機動目標的高階狀態(tài)模型有助于提高估計精度。根據(jù)目標狀態(tài)建模思想的不同，方法大體可分為2類：①對目標未知加速度或加加速度建模[3-4]，構(gòu)建CA、Singer和Jerk等運動學模型，這種方法通用性較強，具有較快的機動響應能力，但估計精度有限；②對目標未知氣動力建模[5]，通過狀態(tài)增廣[6-9]或輸入估計[10-12]等方法構(gòu)建動力學模型，這種方法物理含義清晰，穩(wěn)定跟蹤時估計精度高，廣泛應用于機動再入目標跟蹤。缺點是當氣動力變化較快或模型不夠準確時，跟蹤精度會大幅下降。由于HGV運動形式復雜，單一運動模型難以準確刻畫其運動特征，交互多模型(Interacting Multiple Model, IMM)方法可以較好地表現(xiàn)HGV的復雜機動[8-9,12-13]。然而，在面對狀態(tài)維度大、傳感器采樣率高的HGV估計問題時，IMM往往會大幅增加了運算負擔。

據(jù)以上分析可知，對未知氣動力準確建模是提高HGV跟蹤精度的有效途徑。同時，IMM方法可以顯著改善模型匹配程度，但如何調(diào)和模型集覆蓋能力與計算量之間的矛盾是亟需解決的問題。為此，本文首先論證氣動加速度建模的必要性，將氣動加速度在轉(zhuǎn)彎和俯仰方向存在的耦合信息引入氣動參數(shù)模型，構(gòu)建滾轉(zhuǎn)(Bank)和螺旋(Sprial)模型。然后，綜合考慮估計精度和實時性的需求，使用一種目標狀態(tài)與氣動參數(shù)分離估計的模型替代傳統(tǒng)的狀態(tài)增廣模型，并對大幅變化的參數(shù)對應的機動頻率進行IMM建模。最后，對不同條件下的多種模型跟蹤算法進行仿真驗證。

1 目標氣動加速度分析

(1)

式中：aA為氣動加速度，是造成目標機動的原因；aG為重力加速度。

從式(1)中可以看出，目標在重力加速度aG作用下，運動方式表現(xiàn)為慣性彈道再入運動；在aA各向分量的作用下，目標具有非慣性的跳躍再入和橫向機動能力。為從動力學角度理解aA的變化規(guī)律，將aA在位于目標質(zhì)心的半速度(Velocity-Turn-Climb, VTC)坐標系中分解[5]：

aA=-Duv+L(-utsinφ+uccosφ)

(2)

式中：uv、ut和uc為VTC坐標系中各坐標軸的單位矢量;φ為滾轉(zhuǎn)角；D和L分別為阻力加速度和升力加速度，其表達式為

其中：m為質(zhì)量；ρ為大氣密度；v和S分別為目標速度和等效截面積；α為迎角，φ和α為HGV的控制變量；CD(α)和CL(α)分別為阻力系數(shù)和升力系數(shù)，為α的函數(shù)。

從式(2)可以看出，HGV通過調(diào)整φ和α，使aA各向分量發(fā)生變化從而實現(xiàn)各種復雜機動。為進一步分析aA對目標運動狀態(tài)的影響，對aA進行求導，則有[14]

(3)

2 耦合氣動參數(shù)建模

2.1 氣動參數(shù)定義

HGV機動是由目標的控制變量變化決定的。據(jù)式(2)可知，即使φ和α不發(fā)生變化，由于目標跳躍飛行的原因，大氣密度ρ劇烈變化，氣動加速度aA的各向分量會大幅波動。因此，直接對氣動加速度建模會造成aA的變化規(guī)律難以刻畫。利用已知的狀態(tài)依賴參量ρv2/2將氣動加速度aA中的未知參數(shù)(CD(α)、CL(α)、S、m和φ)轉(zhuǎn)化為氣動參數(shù)αP=[αv,αt,αc]T是一種流行的建模做法[5]。參照文獻[5]建立氣動加速度aA與氣動參數(shù)αP的關(guān)系式為

(4)

氣動參數(shù)αP與大氣密度ρ變化無關(guān)，可以直接表征目標的各向機動：αt>0和αt<0分別對應目標水平左轉(zhuǎn)和右轉(zhuǎn)機動；αc>0和αc<0分別對應目標爬升和俯沖機動。利用氣動參數(shù)α取代氣動加速度aA表示目標的未知機動輸入，可以直接從目標機動特點出發(fā)調(diào)整模型參數(shù)，有助于提高估計精度。然而，αP=[αv,αt,αc]T本質(zhì)上是對氣動加速度aA的解耦處理，分析式(2)可知，aA在轉(zhuǎn)彎和俯仰方向存在明顯的耦合關(guān)系。為有效利用這一耦合關(guān)系以提高建模精度，2.2節(jié)和2.3節(jié)將其作為先驗信息引入氣動參數(shù)αP表達式中，分別建立Bank和Spiral模型。

2.2 Bank模型

聯(lián)立式(2)和式(4)，可得氣動參數(shù)αP與控制變量α、φ之間的關(guān)系式為

(5)

從式(5)中不難看出滾轉(zhuǎn)角φ通過改變CL(α)S/m在縱向平面的指向?qū)崿F(xiàn)目標機動，表現(xiàn)為轉(zhuǎn)彎力參數(shù)αt和爬升力參數(shù)αc的非線性負相關(guān)：

(6)

(7)

氣動參數(shù)αB從動力學的角度表征了目標的控制規(guī)律，具有清晰的物理意義：αv和α1對應阻力和升力的大小，φ表示升力的方向。通過這種形式的參數(shù)建模，αB考慮了αt和αc之間的耦合關(guān)系，但本質(zhì)上與αP同階，并未增加計算復雜度。聯(lián)立式(4)和式(7)，可利用αB表示氣動加速度aA為

(8)

采用一階馬爾可夫模型表示各參數(shù)的機動特性，氣動參數(shù)αB的系統(tǒng)模型可表達為

(9)

式中：wv、wl和wφ分別為對應參數(shù)的過程噪聲；λv、λl和λφ分別為對應參數(shù)的機動頻率。HGV飛行過程中，α主要起到維持目標穩(wěn)定飛行的作用，α頻繁調(diào)整容易造成目標氣動熱急劇增大甚至飛行失穩(wěn)，一般取為常值或小區(qū)間變化[15]。同時，研究表明短時間內(nèi)CD(α)和CL(α)隨α呈現(xiàn)近似線性變化[16-17]。因此，λv和λl可取為較小的值。同時，作為目標機動的主要控制變量，滾轉(zhuǎn)角φ決定了升力在轉(zhuǎn)彎和爬升方向分量的大小，變化較為頻繁，λφ可取為較大的值。

2.3 Spiral模型

與Bank模型思想類似，Spiral模型也考慮了αt和αc之間的耦合關(guān)系。根據(jù)式(5)分別對αt和αc求導，則有

(10)

可知αt和αc可組成一個穩(wěn)定的諧波振蕩器[6]

(11)

式中:?為滾轉(zhuǎn)角φ的角加速度，反映了滾轉(zhuǎn)角φ變化的劇烈程度。將?增廣到氣動參數(shù)αP中，可得改進的氣動參數(shù)向量為αS=[αv,αt,αc,?]T。不難看出，αS利用αv、αt和αc表征了目標在VTC坐標系中各向的機動，利用?表征了αt和αc之間的耦合關(guān)系，?本質(zhì)上是αP的一階導數(shù)信息。

Spiral模型的氣動加速度aA表達式與式(4)保持一致。采用一階馬爾可夫模型表示?的機動特性，氣動參數(shù)αS的系統(tǒng)模型可表達為

(12)

式中：wt、wc和w?分別為對應參數(shù)的過程噪聲；λ?為?對應的機動頻率。由于?變化是目標機動的主要影響，λ?可取為較大的值。

3 IMM分離估計算法

3.1 分離估計模型

基于氣動參數(shù)建模的跟蹤算法通常將氣動參數(shù)增廣到狀態(tài)模型進行濾波計算，但這種算法在多模型條件下會顯著增加計算復雜度?？紤]到需要對氣動參數(shù)進行多模型假設(shè)，本文使用一種狀態(tài)量與氣動參數(shù)分離估計模型替代增廣模型，可有效降低計算復雜度[18]。

圖1 分離估計模型結(jié)構(gòu)Fig.1 Structure of decomposed estimation model

3.2 濾波器描述

3.2.1 狀態(tài)濾波器

ENU坐標系下HGV狀態(tài)方程的離散化表達式為

x(k+1)=Fxx(k)+

G[aA(k)+aG(k)]+wx(k)

(13)

式中：

(14)

(15)

其中：Re為地球半徑。

在球坐標系中建立量測方程，目標的量測值為z(k)=[R(k),A(k),E(k)]T，R(k)、A(k)和E(k)分別為觀測的徑向距離、俯仰角和方位角。根據(jù)ENU坐標系和球坐標系的轉(zhuǎn)換關(guān)系，目標的量測方程可表達為

z(k+1)=h[x(k+1)]+vx(k+1)=

(16)

式中:vx為狀態(tài)濾波器的量測噪聲;nR(k+1)、nA(k+1)和nE(k+1)為各量測值對應的量測噪聲。狀態(tài)濾波器的過程噪聲Qx(k)與量測噪聲Rx(k)參考文獻[19]。

3.2.2 Bank參數(shù)濾波器

對式(9)離散化，可得Bank參數(shù)濾波器狀態(tài)模型的表達式為

(17)

(18)

(19)

(20)

利用狀態(tài)濾波器得到的狀態(tài)誤差協(xié)方差Px(k)求解量測誤差協(xié)方差RB(k+1)為

RB(k+1)=diag(rB1,rB2,rB3)

(21)

式中:

3.2.3 Spiral參數(shù)濾波器

對式(12)離散化，可得Spiral參數(shù)濾波器狀態(tài)模型表達式為

(22)

(23)

vS(k+1)

(24)

根據(jù)狀態(tài)濾波器得到的狀態(tài)誤差協(xié)方差Px(k)求解量測誤差協(xié)方差Rs(k+1)為

Rs(k+1)=diag(rS1,rS2,rS3,rS4)

(25)

式中:

式中：σ?為?對應的量測誤差協(xié)防差。

3.3 IMM設(shè)計

基于氣動參數(shù)建模的跟蹤算法的共性缺點是大多都基于某種單一的假設(shè)模型，一旦在跟蹤過程中目標軌跡和設(shè)定的模型不匹配，或者氣動參數(shù)變化導致加速度狀態(tài)突變時，狀態(tài)估計精度就會迅速惡化。當前普遍認為IMM是解決此類問題的有效途徑。

在氣動參數(shù)系統(tǒng)模型式(17)和式(22)中，氣動參數(shù)時間序列相關(guān)性的表現(xiàn)為機動頻率λ的大小，其取值決定了HGV機動時的狀態(tài)估計精度[8]。當目標穩(wěn)定飛行時，量級為10-4s-1的λ可表現(xiàn)參數(shù)的緩變特性；對于目標機動飛行的情況，氣動參數(shù)瞬時變化較大，量級為0.1～1 s-1的λ可表現(xiàn)參數(shù)的突變特性。同時，在飛行過程中HGV的機動模式可能介于緩變和突變之間任意模式。為此，可利用多個機動頻率λ組成的IMM模型提高模型匹配程度。為調(diào)和機動模式的覆蓋需求與計算復雜度的矛盾，本文結(jié)合上述取值范圍，對Bank模型的φ和Spiral模型的?對應的機動頻率λφ和λ?進行適當選取。如圖2所示，設(shè)計的IMM模型由設(shè)置r個不同機動頻率λi(i=1,2,…,r)的濾波器組成。

圖2 氣動參數(shù)IMM示意圖Fig.2 IMM schematic of aerodynamic parameter

4 仿真分析

為驗證本文方法可行性，設(shè)計如下仿真環(huán)境：①目標參數(shù)：HGV模型參考美國洛馬公司CAV-H的基本參數(shù)[20]，飛行器總長約為4 m，質(zhì)量1.7 t，初始高度45 km，速度為4 000 m/s，采用三自由度動力學方程積分生成仿真彈道，采用標準制導法控制φ和α變化實現(xiàn)機動[15]；②傳感器參數(shù)：采用擴展卡爾曼濾波(EKF)濾波器，距離量測標準差為500 m，方位角、俯仰角量測標準差均為0.01 rad，Bank模型和Spiral模型中r取值均為3，對應的機動頻率取值分別為0.001、0.01和0.5 s-1。為檢驗算法性能，對仿真算例進行500次蒙特卡羅仿真，分別計算均方根誤差(Root Mean Square Error, RMSE)。

4.1 不同模型條件下估計性能仿真

利用2類模型跟蹤算法進行仿真對比：①運動學模型：CA、Singer和Jerk模型；②動力學模型：基于氣動參數(shù)αP=[αv,αt,αc]T增廣的機動載入飛行器(MaRV)模型、本文Bank模型和Spiral模型。需要注意的是，為保證仿真合理性，仿真中的6個模型均為單模型，即不考慮IMM情況。

圖3 不同模型位置和速度估計RMSEFig.3 RMSE in position and velocity estimate of different models

圖3給出6種模型位置與速度RMSE仿真結(jié)果。不難看出：①當目標發(fā)生機動時(氣動參數(shù)變化)，各模型的位置和速度跟蹤誤差發(fā)生顯著增大，然后逐漸減??；②動力學模型跟蹤誤差收斂較快，顯著小于運動學模型，表現(xiàn)出較好的跟蹤精度；③在動力學模型中，Bank模型性能優(yōu)于其他2種模型，Spiral模型略微優(yōu)于MaRV模型，但差別不大。

4.2 不同動力學模型條件下IMM估計性能仿真

為分析氣動參數(shù)估計性能，對不同動力學模型條件下IMM算法進行仿真，MaRV模型的機動頻率參考第4節(jié)仿真環(huán)境數(shù)據(jù)。圖4和圖5分別給出3種動力學模型氣動參數(shù)估計值以及相應的RMSE仿真結(jié)果。其中，Bank模型的αt和αc的估計值是利用式(7)間接計算得到。從圖4和圖5中可以看出：①在當氣動參數(shù)劇烈變化時，不同動力學模型對氣動參數(shù)的估計均發(fā)生時延，估計誤差逐漸增大，且難以快速收斂，說明動力學模型難以精確估計氣動參數(shù)；②MaRV模型與Spiral模型對氣動參數(shù)的估計誤差量級相當，而Bank模型的估計精度則優(yōu)于上述2種動力學模型。

圖4 不同動力學模型氣動參數(shù)估計值Fig.4 Aerodynamic parameter estimation values of different dynamic models

圖5 不同動力學模型氣動參數(shù)估計RMSEFig.5 RMSE in aerodynamic parameter estimate of different dynamic models

對動力學模型仿真結(jié)果的原因進行分析：①MaRV模型未考慮各向氣動參數(shù)耦合的先驗信息，是一種相對“寬松”的解耦模型，導致估計性能下降；②Spiral模型中?本質(zhì)上是加速度的一階導數(shù)信息，狀態(tài)濾波器無法提供關(guān)于?的先驗信息，其僅能通過Spiral參數(shù)濾波器間接辨識，導致?估計精度有限；③Bank模型中αB本質(zhì)上加速度同階，因此利用狀態(tài)濾波器計算得到αB的偽量測值，由于考慮了氣動參數(shù)耦合信息，Bank模型能夠有效提高估計精度。

4.3 不同濾波器參數(shù)條件下IMM估計性能仿真

為檢驗不同濾波器參數(shù)對本文所提動力學跟蹤算法的影響，利用表1中6個算例從濾波算法、噪聲方差和采樣周期對模型性能進行分析。

利用算例1和2對不同濾波算法下動力學模型進行性能仿真，如圖6所示為位置與速度RMSE仿真結(jié)果。不難看出，無跡卡爾曼濾波(UKF)的估計精度雖然高于EKF，但相對于動力學模型對估計精度的提升，濾波算法的優(yōu)勢并不明顯。可見，模型的精細程度對跟蹤算法影響更大。

利用算例1和算例3～6對不同量測參數(shù)條件下的動力學模型進行性能仿真，如表2所示為位置、速度和氣動參數(shù)的RMSE均值對比情況?？梢钥闯觯孩匐S著量測噪聲方差增大或減小，各狀態(tài)量RMSE均值均隨之增大，但氣動參數(shù)估計誤差變化更為劇烈，可見量測噪聲對氣動參數(shù)估計影響較大；②隨著采樣周期增大或減小，各狀態(tài)量RMSE均值隨之增大或減小，但位置與速度估計誤差變化更為劇烈，可見采樣周期對位置與速度估計影響較大。

分析以上仿真結(jié)果可以得出以下結(jié)論：①濾波算法對不同動力學模型跟蹤算法性能影響有限，采用EKF可有效提高計算效率；②傳感器的噪聲方差大小直接影響了算法對目標狀態(tài)的估計精度；③調(diào)和氣動參數(shù)辨識精度與計算復雜度的矛盾，應當合理選擇采樣周期；④雖然濾波算法和和量測參數(shù)的變化使模型性能發(fā)生了改變，但總體上Bank性能優(yōu)于Spiral模型，Spiral模型略微優(yōu)于MaRV模型，進一步證實了4.1節(jié)和4.2節(jié)的結(jié)論。

表1 性能對比算例參數(shù)設(shè)置

算 例位置RMSE均值/m速度RMSE均值/ (m·s-1)αvRMSE均值/10-6αtRMSE均值/10-6αcRMSE均值/10-6MaRVSpiralBankMaRVSpiralBankMaRVSpiralBankMaRVSpiralBankMaRVSpiralBank1107.8105.797.0023.9423.0019.347.0696.7701.86180.8577.9571.2219.8620.0617.503173.2165.6158.532.8729.5425.8910.667.7892.50084.2281.2874.3020.5821.6418.114192.0187.4171.337.1232.6526.809.1328.5852.04786.5183.7876.3721.8222.1818.35585.6680.6176.6020.0119.4717.964.2274.1061.34374.2173.2668.3018.5819.6417.206167.6160.3142.633.5932.8926.4818.0417.4912.56100.797.7089.9723.2523.3219.41

5 結(jié) 論

1) 對HGV的加速度導數(shù)項進行了推導和分析，得出機動的實質(zhì)是氣動加速度改變引起目標質(zhì)心運動和繞質(zhì)心運動變化的結(jié)論。為改善算法估計精度，應對目標的高階狀態(tài)進行建模。

2) 在非耦合氣動參數(shù)建模的基礎(chǔ)上，將氣動參數(shù)在轉(zhuǎn)彎和俯仰方向存在的耦合關(guān)系作為先驗信息，構(gòu)建了Bank和Spiral模型。

3) 針對多模型條件下增廣模型導致計算復雜度增加的問題，使用一種由狀態(tài)濾波器和氣動參數(shù)濾波器組成的分離估計模型替代增廣模型。同時，對模型中變化較大的參數(shù)進行機動頻率IMM建模，提高了模型匹配程度。

4) 仿真結(jié)果表明，動力學模型性能顯著優(yōu)于運動學模型。同時，在不同條件下對比3種動力學模型可知，本文所提的2種模型對位置、速度和氣動參數(shù)的估計精度好于MaRV模型，Bank模型顯著優(yōu)于Spiral模型，且計算復雜度更小。下一步將繼續(xù)深入研究HGV動力學特性，以期構(gòu)建更加精細的動力學模型，進一步改善估計精度。