臨近空間高超聲速目標跟蹤動力學模型

李 凡，熊家軍，張 凱，韓春耀，席秋實，張懷念

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

0 引 言

臨近空間高超聲速飛行器(Near space hypersonic vehicle, NSHV)一般指飛行于20～100 km高空、飛行馬赫數大于5的飛行器[1]，隨著NSHV在軍事領域發展及應用，其巨大的戰略價值受到全世界的廣泛關注，西方各軍事大國投入大量資源進行NSHV關鍵技術的研究。2017年10月30日美國于夏威夷考艾島導彈試驗基地成功進行了潛射型高超聲速助推滑翔導彈的首飛試驗(試驗代號CPS FE1)，此后，美國海軍戰略系統項目辦公室(SSP)聲稱，已著手在上裝載這一新型的垂發模塊，“弗吉尼亞”級攻擊型核潛艇將成為這一常規快速打擊武器的發射平臺，這一舉措預示著美國NSHV已初步具備實戰能力，有分析認為美國有望在未來幾年啟動NSHV裝備采辦程序。NSHV快速發展提高了進攻方的軟殺傷及硬摧毀能力，但對于防御方而言，各國現有的預警監視體系對這類飛行器跟蹤問題還沒有很好的解決方案，因此，進行NSHV目標跟蹤的關鍵技術研究具有重要意義。

機動目標跟蹤問題實質上是對未知系統的描述，其單目標跟蹤中主要包括濾波算法與機動模型兩個方面[2]。濾波算法聚焦于局部量測噪聲的誤差修正，其本質是通過多個時刻殘差樣本提取量測噪聲的分布特征參數，自適應修正當前時刻目標的狀態估計值，濾波算法受噪聲環境影響較大，對不同運動模式具有較強的適應性，此外，濾波算法發展依賴于基礎數學理論的研究，難度相對較大；當前跟蹤算法研究集中于機動模型改進或創新，針對不同機動樣式目標跟蹤的區別主要在于機動模型建立，機動模型致力于目標運動特性的描述，通過建立合理模型逼近目標真實運動模式，在觀測數據有限的情況下，好的模型對于機動目標跟蹤問題至關重要。當前高速強機動目標機動建模方法大體分為兩條主線，一是從運動幾何的角度研究目標狀態隨時間變化的規律[3-4]，直接以加速度(或加加速度)為基點，分析加速度之間的相關特性，并建立線性離散化的運動學模型跟蹤狀態方程，通過對加速度在時間維的積分求和得到目標狀態估計，因此建模過程中不涉及坐標變換，該方法將目標看作一個等效質點，不考慮目標加速度產生的原因、目標與地球之間的引力關系及目標的姿態變化。常用的模型包括單模型及多模型，現有單模型主要有勻速(Constant velocity, CV)、勻加速(Constant acceleration, CA)、Singer、“當前”統計模型(Current statistical, CS)、Jerk[5-6]，這類模型結構簡單，計算量小，但對復雜運動及強機動運動適應性較差；多模型主要包括固定結構多模型(Fixed structure multiple model, FSMM)和變結構多模型(Variable structure multiple model, VSMM)[7-8]，模型集覆蓋的運動模式廣，機動適應性較強，但噪聲模型選取及模型交互參數設置較為困難。總體而言，運動學模型狀態量之間具有很好的線性關系，通過狀態轉移矩陣遞歸實現狀態預測，先驗信息需求較弱，這一方法實質是對目標運動模式的趨近。二是從目標受力的角度出發[9]，對目標進行力學分析推導各方向的加速度特性，該方法分別在不同的坐標系下考慮重力、空氣動力及推力(有動力的情況)的影響，由于目標受力分析分別在不同的坐標系下表達較為簡潔，需要轉換到同一坐標系下進行目標狀態量估計及濾波，因此模型的非線性程度較高，此外，考慮力產生機動大多存在較為嚴重的耦合，該方法優點在于從彈道設計角度考慮目標機動的產生，模型建立較為合理時能較為準確地估計氣動參數，但先驗信息要求較強。現有NSHV目標跟蹤動力學模型存在兩個方面的問題，其一是動力學模型大多只考慮了空氣動力因素即無動力情況，但NSHV存在有動力巡航方式，該方式飛行器控制性能、機動性能更佳，同時跟蹤難度也更大。其二是動力學模型大多假設氣動力參數服從為白噪聲或一階馬爾科夫過程，需要建立更加合理的模型描述其氣動加速度特性。

本文針對跳躍式NSHV目標提出了一種新型的動力學跟蹤模型，從目標運動特性出發，對氣動參數建立了新模型，使之同時體現周期性與衰減性，同時，為增強算法對機動突變點的適應性，引入強跟蹤濾波進行加速度突變檢測及濾波修正。

1 NSHV運動特性分析

平衡滑翔及跳躍式滑翔是NSHV最為常用的兩種飛行方式，其中跳躍式滑翔能夠在縱向平面進行多次跳躍，機動性能較強，是當前機動目標跟蹤領域的一大難題，本文僅討論跳躍式NSHV的跟蹤。而對目標運動特性進行分析并抽象出其運動模式的一般規律，是構建合理跟蹤模型的前提條件，因此，本節主要分析了NSHV在不同維度的運動特點以及可能的機動樣式。

1.1 NSHV動力學分析

飛行器動力學分析有利于理解NSHV跳躍式飛行原理、分析其運動特性，當飛行器在空間內飛行時一般都有三個力作用于該物體上，即重力、推力、空氣動力，即飛行器受力表達式為[10]：

F=FT+FA+FG=FT+FD+FL+FG

(1)

其中,FT為發動機推力(無動力時推力為0)，FA為空氣動力主要分為空氣阻力FD及升力FL，FG為重力。則NSHV在助推段的離散化動力學方程如下：

式中：k為時間刻度，ν(k)為目標速度，α(k)為攻角，M(k)為機體質量，g為重力加速度，θ(k)為航跡傾角，σ(k)為航向角，r(k)為目標質心到地心距離，φ(k)為目標緯度。

在助推段分析推力及空氣動力對加速度dν(k)/dk及縱向平面角速度dθ(k)/dk的影響程度，分別假定式(2a)及式(2b)中除推力及空氣動力外參數為常數，易知在FT(k)cosα(k)與空氣阻力FD(k)在加速度dν(k)/dk中的靈敏性相同，同理，FT(k)sinα(k)與空氣升力FL(k)在dθ(k)/dk中的靈敏度一致，此外，推力在加速度及縱向平面角速度中分別為正弦和余弦對偶分量，因此，推力FT(k)至少有一個分量與空氣動力分量在助推段的靈敏性相同。

NSHV目標助推拉升后進入無動力滑翔階段，這一方式能有效規避雷達探測，降低視距。其中無動力滑翔階段動力學方程如下：

(3)

式中：φ(k)為傾側角，當目標無動力飛行時受力較為簡單，只存在空氣動力與重力，飛行器通過自身調整實現橫向機動。目標航向角σ(k)較小[11](σ(k)∈[0°,10°])，即目標在橫向機動強度較小；根據文獻[11]可知航跡傾角變化率dθ(k)/dk較大，目標在高度層機動強度較大。

當有動力與無動力互相轉換時相當于加入一個階躍的推力，根據分析可知至少有一個推力分量與空氣動力對運動影響程度處于同一量級，僅考慮無動力情況時，動力學模型難以描述轉換階段的加速度突變，轉換階段容易出現跟蹤誤差升高。

1.2 橫縱向機動特性分析

縱向平面內“打水漂”跳躍飛行是這類目標最鮮明的特點，其軌跡跳躍幅度呈衰減形式，類似于振蕩運動，可以推測相鄰整數周期的目標狀態量呈現類周期性變化[12-15]。

從飛行器控制參量方面看，攻角與傾側角控制飛行器的主控量，其中攻角控制機體縱軸與速度方向間的俯仰幅度，通過不同的升力/阻力系數決定目標在縱向上的飛行高度、飛行速度、航跡傾角等，攻角的設定需要滿足飛行走廊的限制條件，不可能非常復雜，攻角經常采用常值、分段或線性等簡單函數來描述。傾側角同時影響目標在橫縱向的飛行軌跡，決定了飛行器的橫向機動，為了有效規避、突防禁飛區，傾側角設定通常是在飛行走廊內幅度不變符號反轉。

攻角方面。文獻[16]指出在滿足飛行走廊條件限制前提下無論如何設定攻角控制量，只要不是一直滿足平衡滑翔條件，其彈道都呈現出跳躍滑翔樣式。實際上平衡滑翔比跳躍滑翔更難實現，平衡滑翔方式對控制參量設定要求更高，一般來說，跳躍滑翔攻角取值可以是滿足飛行走廊條件下的一個區間，而平衡滑翔攻角取值只能是滿足式(4)條件的一個定值。

(4)

平衡滑翔狀態下，航跡傾角很小，一般可以認為cosθ(k)=1。其中

FL(k)=CL(k)q(k)S

(5)

q(k)=ρ(k)v2(k)/2=ρ0e-βh(k)v2(k)/2

(6)

式中:CL(k)為升力系數，是與控制量攻角α(k)和馬赫數Ma相關的函數，q(k)為動壓，S為飛行器的參考面積，ρ(k)為當地大氣密度，ρ0為海平面大氣密度，β為高度常數，h(k)為目標高度。

文獻[16]指出當其他參數一致僅改變初始高度時，若飛行器滑翔初始高度不滿足平衡滑翔狀態時，NSHV軌跡都會產生跳躍。當初始高度h2(k)小于平衡滑翔高度h1(k)時，當地大氣密度ρ0e-βh2(k)大于平衡滑翔的ρ0e-βh1(k)，導致飛行器升力增加，無法滿足平衡滑翔條件，其合力作用向上，NSHV向上做加速運動，當地大氣密度降低，升力減小直到合力為0，向上的速度達到最大，然后向上做減速運動，爬升到最高點后，在重力的作用下向下做加速運動，高度降低，升力增加并再次達到平衡時，向下的速度最大，高度繼續降低直至向下速度為0，飛行器到達最低點，此后繼續向上飛行，完成一個軌跡周期。(初始高度大于平衡滑翔高度也會形成跳躍滑翔軌跡)。此外，僅改變初始速度或初始速度傾角時，其彈道同樣為跳躍滑翔樣式。因此，NSHV升力縱向的升力加速度相關性也具有類似的周期特性。

傾側角方面。傾側角的設定主要是規避敵方探測攔截，美國最初將NSHV目標定義為全球快速打擊武器，旨在高速對敵方大縱深目標進行軟殺傷或硬摧毀，可以通過大C或S形橫向彈道規避敵方重點防御區域后依靠高速打擊縱深目標，橫向高頻跳躍或其他方式復雜機動會導致在敵方區域的無效路徑增長，被發現的概率增大，攻擊的時效性降低。此外，NSHV飛行速度較快，橫向機動半徑較大，飛行

器動能有限，無法如縱向進行多次跳躍。

2 動力學跟蹤模型構建

假定雷達站坐標為東-北-天坐標系(East-north-up, ENU)，在半速度(Velocity-turn-climb, VTC)中分析氣動力模型及推力模型，并在地心慣性坐標系(Earth central inertial, ECI)中構建重力及表視力模型[10]。

(7)

2.1 重力及表視力模型

假設目標在ECI及ENU坐標系下位置狀態分別為XECI,XENU，假定雷達站的經緯高為L,B,H，可知ECI到ENU的轉換為[5]：

(8)

其中轉換矩陣為

Re為地球等效半徑，ωet為從參考時刻到t時刻ECI和地心坐標系(Earth central, EC)的夾角(一般假定雷達探測到目標開始EC和ECI坐標系重合)[17]，則ENU坐標系下目標加速度為：

(9)

2.2 氣動力模型

推力產生加速度是通過對空氣作用而實現的，因此，將推力模型認為是一種疊加的特殊空氣動力，則aP=aT+aA可表示為：

(10)

根據VTC坐標系與ENU坐標系的轉換規則，ENU坐標系下新的氣動加速度為[17]：

(11)

2.3 氣動加速度建模

目標在飛行過程中重力加速度變化較小，其運動模式變化主要取決于氣動加速度的改變，第2.2節中分析了各項加速度的物理含義，結合第1.2節中NSHV機動特性分析，對阻力加速度，橫向轉彎加速度及爬升加速度進行分布假設建模。

1)對于爬升加速度，由縱向機動特性可知，其自相關存在周期性與衰減性，因此，采用阻尼函數表示爬升加速度的相關特性。其自相關函數為：

(12)

1/[α+(β-w)j]+1/[α-(β+w)j]}

(13)

其中，F(·)為傅里葉變換，w為角速率，令G(w)=[α-(β-w)j][α+(β+w)j]，則有

W(wj)H(wj)H(-wj)

(14)

其中，W(wj)為白噪聲的傅里葉變換。則濾波器的傳遞函數為：

(15)

其中，s為拉普拉斯算子，則機動加速度a(t)的連續時間微分方程為：

(16)

基于動力學模型跟蹤NSHV，需要將VTC坐標系下的氣動加速度進行轉換到ENU坐標系下，這一過程存在高度的轉換耦合，因此難以得到過程噪聲協方差的顯性表達式，但狀態協方差的迭代更新可以一定程度適應過程噪聲協方差取值變化，所以更為關注狀態轉移矩陣的離散化，式(16)可簡化為：

(17)

2)對于轉彎力加速度，根據橫向機動特性分析可知，其機動強度相對縱向較弱，常用的彈道形狀為直線、C形或S形，其橫向加速度相關性呈較強的衰減性，因此假設橫向加速度服從一階馬爾科夫模型，轉彎加速度的微分方程為：

(18)

式中：γ為機動時間常數。

3)對于阻力加速度，當目標質量、等效面積等恒定時，目標所受阻力僅與速度、空氣密度有關。

(19)

根據此前分析可知，跳躍滑翔式目標在縱向上的高度呈類周期的跳躍形式，而目標總體速度呈階梯下降的趨勢，因此，其阻力加速度也呈現振蕩下降的形式。因此，將其建模為衰減振蕩自相關隨機過程。

3 動力學跟蹤模型狀態方程

為方便討論，在此區分運動狀態量與氣動狀態量的概念，運動狀態量為ENU坐標系下狀態量(包括位置、速度、加速度)，氣動狀態量為VTC坐標系下狀態量(包括加速度及加加速度)。

3.1 運動狀態量更新方程

k時刻運動狀態量更新方程fI(X(k))為[17]：

(20)

3.2 氣動狀態量更新方程

式(18)的離散化方程為：

at(k+1)=e-γTat(k)+wt(k+1)

(21)

其中，T為采樣時間間隔。式(17)的離散狀態方程可通過式(22)求得。

(22)

(23)

式中:

則完整的跟蹤狀態離散方程為：

4 加速度突變點的濾波修正

根據式(2)～(3)可知，有無動力段轉換時目標產生較大的機動加速度，將導致跟蹤誤差增大，在跟蹤初始模型經過一段時間的調整與目標運動模式匹配，之后較為穩定地描述目標運動，并通過增益進行預測量的修正，模型需要保證跟蹤算法整體的魯棒性，其加速度不會出現突變，當目標出現機動突變時模型無法進行瞬態變化，同時，濾波算法修正誤差需要一定的量測樣本，因此會出現機動跳變位置的誤差升高。針對這一問題引入強跟蹤濾波對機動突變點進行補償修正。考慮到有動力轉換主要用于推動飛行器從彈道低點拉升，選用升力加速度ac(k)作為機動突變檢測參數。判決量：

(24)

當相鄰加速度相差達到兩個數量級以上即ξ(k)≥2時，認為發生機動突變，采用強跟蹤濾波STF對當前時刻量測進行濾波修正[18]。

5 仿真分析

圖1 仿真軌跡

5.1 仿真試驗1

航跡I為有動力巡航軌跡，航跡I初始位于助推段從雷達視距中出現，達到助推頂點后經過一個完整的滑翔—助推，完成跳躍巡航。三種算法的跟蹤仿真結果如圖2所示。

圖2 位置及速度誤差均方根比較

如圖2所示，算法I直接通過運動加速度描述目標的位置及速度變化，在濾波初始階段,目標運動加速度變化的量級遠超過氣動加速度，算法I需要一定的時間適應加速度變化速度，其位置誤差及速度誤差收斂速度較慢；在200～300 s內由于運動加速度出現突變，導致位置及速度濾波誤差都有一定程度躍起，其機動適應性相對較差。其余時刻目標機動較弱，運動模式趨于勻速及勻加速之間，三種算法位置濾波精度相近，但該算法速度濾波精度較高，略優于算法II，與算法III相當。

算法II基于氣動加速度擴維實現目標狀態預測及濾波，氣動加速度變化的量級小，但本身的量級也較小(參見圖3中氣動加速度量級)，因此，需要一段時間才能穩定描述氣動加速度變化，這是限制動力學模型收斂時間的關鍵性因素之一；在200～300 s內強機動時其位置及速度誤差相對算法I較小，動力學模型對運動加速度的突變具有一定的適應性；其余時刻起誤差水平較為平緩。

算法III基于衰減振蕩函數構建新型動力學跟蹤模型，并通過強跟蹤濾波檢測氣動加速度突變點實時修正濾波狀態，其收斂速度較快；在200～300 s內速度濾波誤差起伏較小，但位置誤差保持平穩趨勢。整體而言本文算法整體水平較為穩定，濾波效果最好。為更好地比較算法性能，給出各算法氣動加速度濾波估計仿真圖。(其中算法I的氣動加速度估計是通過氣動加速度到運動加速度逆運算實現)。

圖3 氣動加速度參數估計

在氣動加速度方面，算法I通過濾波狀態解算得到氣動加速度量，如圖3所示，算法I的氣動加速度估計最大特點為過調整幅度較大，因運動加速度變化量較大導致調整速度較快、過調整較為嚴重，此外，在200～300 s內，由于運動加速度在這一時間段內變化較大，該算法解算的氣動加速度在此時間段內出現局部的大幅度跳躍，但該時間段內目標實際氣動加速度并沒有劇烈變化。

算法II氣動加速度量的估計明顯優于算法I，過調整量相對較小，但調整時延較長，在200～300 s內，能較為準確反應真實氣動加速度的變化趨勢。

算法III采用了新型模型，對氣動加速度量估計效果最好，其包絡與真實值匹配程度最高，此外，采用氣動加速度檢測及強跟蹤濾波補償，其時延相對較小。

圖4 真實運動加速度

如圖3所示，氣動加速度方面在0～100 s內快速變化，其余時刻曲線較為平緩；而運動加速度方面如圖4所示，出現兩次變化較大，造成兩種加速度不對應的原因分析如下：

1)氣動加速度是運動加速度的分量，由式(20)可知，氣動加速度為運動加速度分量，運動加速度還包括重力及表視力，特別是重力不可忽略。

2)氣動加速度的耦合，由式(11)轉換矩陣可知，氣動加速度分量對運動加速度的影響與當前各方向的速度有關，氣動加速度與運動加速度的轉換存在嚴重的耦合，從圖3及圖4也可以看出,氣動加速度與運動加速度值不在一個數量級，其對應的變化規律難以定性分析。

3)機動的坐標屬性，目標運動都是具有坐標屬性的，在不同的坐標系下運動模式會呈現不同的形式，本文仿真的雷達觀測位置也一定程度造成了氣動加速度與運動加速度的不對應關系。

因此，三種算法位置及速度誤差在200～300 s段出現一定的誤差起伏。無論動力學或運動學跟蹤模型最終都必須過渡到運動加速度、速度進行狀態量估計，運動加速度的突變一定會引起速度及位置急劇變化，這也是氣動加速度平緩而出現濾波誤差突變的原因。

5.2 仿真試驗2

航跡II為無動力滑翔軌跡，跟蹤算法的參數設置不變，得到算法I、算法II及算法III的濾波結果如下圖(篇幅所限，僅給出氣動加速度估計結果)。

圖5 氣動加速度估計

無動力滑翔時，氣動加速度調整量較小，三種算法都能大致反應氣動加速度變化，并呈現局部小范圍的波動，在阻力加速度估計中，算法I及算法II都出現較大波動，算法I容易過調整，估計精度一直處于波動狀態；算法II前期誤差較大，但后期估計精度較為穩定；算法III整體水平最為穩定，精度較高。在轉彎力加速度估計中，算法I同樣存在過調整現象；算法II調整速度較慢，試驗較長；算法III效果較好。對于爬升力加速度估計，前期三種算法估計誤差較大，后期誤差趨于穩定，且三種算法精度相近。

為定量比較算法的濾波性能，位置、速度、氣動加速度誤差統計平均如表1所示。

綜上所述兩個部分的仿真結果，算法I適應勻速到勻加速之間的運動，直接通過反應運動加速度變化，導致對氣動加速度估計容易出現過調整。算法II整體濾波誤差較為穩定，氣動加速度估計略優，但由于氣動加速度量級較小，導致位置誤差收斂速度較慢。將本文算法算法III與算法I及算法II進行比較，可知本文算法所構建模型的優越性。本文主要將氣動加速度建模為衰減振蕩函數形式，對氣動加速度的變化規律匹配程度更高，在位置及速度跟蹤誤差方面精度更高；此外，引入強跟蹤濾波的機動檢測與修正，降低估計時延提高跟蹤精度。

表1 位置、速度及氣動加速度誤差統計平均Tabble 1 Position, velocity and aerodynamic acceleration error statistical average

本文的研究重點在于將氣動加速度自相關假設為具有衰減振蕩的二階馬爾科夫模型，因此，在仿真中僅對比了兩種經典的指數衰減一階馬爾科夫模型，在本文的仿真場景下本文模型跟蹤性能略優，但并不能說明其效果一定優于改進后的經典模型，所幸本文模型從底層的建模假設進行研究，適用于經典模型的改進方法同樣可以運用于本文模型。此外，當前IMM依然是NSHV跟蹤的首選方案，本文模型為后續跟蹤NSHV目標IMM中模型選擇，提供了一個更具潛力的單模型。

6 結 論

本文基于滑躍式NSHV運動方程，結合NSHV控制參量分析了目標的橫縱向機動特性，闡述了滑躍式軌跡產生的原因，在此基礎上構建了一種新型NSHV跟蹤動力學模型，仿真表明本文模型相對現有經典模型有一定優越性，但還存在衰減系數、機動頻率設置，以及在IMM中的應用等其他方面問題有待進一步研究。