張 凱，熊家軍，李 凡，付婷婷

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

0 引 言

作為新一代跨大氣層空天飛行器，以HTV-2、AHW為代表的高超聲速滑翔目標(Hypersonic Gliding Reentry Vehicle, HGRV)結(jié)合航天器與航空器的特征，能利用其高升阻比外形機動形成一個較大的打擊區(qū)域[1-2]。與傳統(tǒng)慣性目標相比，其軌跡預測的可能性大大降低，這給攔截防御帶來巨大挑戰(zhàn)。針對HGRV這類機動不確定目標的軌跡預測問題，采用傳統(tǒng)彈道預測方法難以準確預測，需結(jié)合目標的運動和作戰(zhàn)特點尋求新思路[3]。

現(xiàn)有的HGRV軌跡預測方法大多是從目標歷史飛行數(shù)據(jù)出發(fā)，根據(jù)量測數(shù)據(jù)辨識出目標的氣動力參數(shù)[4-5]、速度傾角[6]、控制變量[7]等參數(shù)，結(jié)合初始估計狀態(tài)積分外推實現(xiàn)軌跡預測。對于防御方而言，這類不考慮目標未來機動的軌跡預測方法能預測的時間窗口十分有限，無法滿足防御作戰(zhàn)需求。雖然HGRV軌跡預測面臨極大困難，但由于目標在高超聲速條件下機動能力受熱防護、氣動布局、制導控制規(guī)律等約束，其機動特征并非完全無章可循[8]。此外，在對空中目標態(tài)勢與威脅評估中，通過推斷目標的飛行意圖可以有效降低未來機動不確定性的影響，這種基于意圖推斷的軌跡預測思路在機器人行為推測[9-10]、民航軌跡預測[11-12]和飛航導彈攔截[13-14]等領(lǐng)域上得到廣泛認可，具有一定借鑒意義。

據(jù)以上分析不難看出，針對HGRV機動不確定問題，結(jié)合目標運動特點合理推斷飛行意圖是提高長期軌跡預測準確性的有效途徑之一。為此，本文首先對HGRV進行動力學建模，利用氣動參數(shù)構(gòu)建HGRV的機動模式集。然后在假定來襲HGRV必定攻擊某目標的前提條件下，結(jié)合其飛行意圖合理構(gòu)造代價函數(shù)，利用貝葉斯原理迭代推導機動模式和運動狀態(tài)的遞推公式，最后通過蒙特卡洛采樣實現(xiàn)貝葉斯軌跡預測算法。期望通過這種基于意圖推斷的貝葉斯軌跡預測思想為HGRV長期軌跡預測提供一定理論指導。

1 問題描述

從圖1中可以看出，狀態(tài)估計反映了目標當前飛行狀態(tài)，推斷的意圖反映了目標飛行目的，僅利用狀態(tài)外推或者飛行意圖對HGRV這類機動再入目標進行軌跡預測都會降低預測精度。本文所提軌跡預測方法旨在結(jié)合兩者優(yōu)勢，這種基于意圖推斷的思路需要解決兩個問題：一是利用動力學模型構(gòu)建HGRV機動模式集，覆蓋目標在未來時刻可能的飛行樣式；二是結(jié)合飛行意圖推斷目標的機動模式，利用機動模式對應(yīng)的預測模型對HGRV進行軌跡預測。下文內(nèi)容分別從這兩點出發(fā)進行闡述。

2 預測模型

機動不確定條件下的HGRV軌跡預測依賴于機動模式集的構(gòu)建。傳統(tǒng)飛航式目標一般利用不同轉(zhuǎn)彎率表示轉(zhuǎn)彎模型表示相應(yīng)的機動模式，由于HGRV的機動能力往往受到特定環(huán)境、技術(shù)因素限制，難以實現(xiàn)類似飛航式目標的復雜機動。因此，HGRV機動模式集應(yīng)當根據(jù)其運動特性設(shè)計，保證預測的軌跡可飛。本節(jié)在動力學建?；A(chǔ)上，結(jié)合目標機動特點，構(gòu)建HGRV的機動模式集。

2.1 動力學模型

未知量氣動加速度a的改變是造成HGRV機動的主要原因。對于HGRV這類機動再入目標的動力學建模問題，比較流行的做法是將氣動加速度a從雷達站東北天(East-North-Up, ENU)坐標系轉(zhuǎn)換到半速度(Velocity-Turn-Climb, VTC)坐標系中，表示成含氣動參數(shù)u=[αv,αt,αc]T的表達式[15]：

(1)

(2)

式中:r為目標瞬時地心距；B為雷達站地理緯度；Re為地球半徑。

2.2 機動模式集

飛行器機動模式由控制規(guī)律決定。HGRV控制參數(shù)包括攻角α和傾側(cè)角υ，攻角α主要作用是維持目標穩(wěn)定飛行，通常取為常值或變化較??；傾側(cè)角υ是實現(xiàn)機動飛行的主要控制參數(shù)，通過改變升力方向?qū)崿F(xiàn)預期的機動動作，表現(xiàn)為爬升力參數(shù)αc和轉(zhuǎn)彎力參數(shù)αt的變化[16]。

已知氣動參數(shù)u與控制變量的關(guān)系式為[15]：

(3)

式中:CD(α)、CL(α)分別為阻力系數(shù)和升力系數(shù)，與攻角α有關(guān)；S為目標等效截面積，m為目標質(zhì)量，均為常值。對于防御方而言，上述參數(shù)與控制變量均為未知量。

再入滑翔過程中，由于攻角α變化較小，且主要影響縱向運動，為減小構(gòu)建機動模式集的復雜度，可忽略攻角α對機動模式的影響，利用不同傾側(cè)角υ對應(yīng)的氣動參數(shù)u表示HGRV在相應(yīng)機動模式下的機動輸入。受到目標氣動特性的限制，υ的變化范圍通?？扇?40～40度，在沒有先驗信息的條件下可假設(shè)υ取值符合正態(tài)分布。為此，在軌跡預測過程中，需要解決氣動參數(shù)u的預測問題。為避免υ變化對氣動參數(shù)u預測的不利影響，定義升力參數(shù)αl：

(4)

結(jié)合式(3)～(4)可知，αv和αl是CD(α)和CL(α)的函數(shù)，與傾側(cè)角υ無關(guān)。研究表明CD(α)和CL(α)隨攻角α近似線性變化[17]，因此，對αv和αl進行預測從而間接得到氣動參數(shù)u的思路是比較合理的。具體步驟：

(1) 參數(shù)估計 根據(jù)量測數(shù)據(jù)利用動力學模型(2)估計出氣動參數(shù)u=[αv,αt,αc]T，通過式(4)計算αl的估計值，并對估計值進行平滑處理。

(2) 參數(shù)預測 可采用以下兩種方法擬合αv和αl變化規(guī)律：①時間序列分析方法：采用(p,d,q)階ARIMA模型[18]；②最小二乘法：構(gòu)造擬合多項式[5]。

(3) 機動建模 根據(jù)不同機動模式下傾側(cè)角υ取值，結(jié)合式(3)～(4)計算相應(yīng)氣動參數(shù)u的預測值。

3 貝葉斯軌跡預測

基于意圖推斷的軌跡預測是根據(jù)實時戰(zhàn)場態(tài)勢進行軌跡規(guī)劃過程：首先是從量測數(shù)據(jù)中提取出對意圖推斷有用的信息，然后按照推理機制對預測軌跡進行規(guī)劃。實質(zhì)是根據(jù)飛行意圖對目標未來狀態(tài)進行估計的問題，即推斷未來某時刻目標在特定空域內(nèi)概率大小的問題。貝葉斯估計能夠有效地利用已經(jīng)獲得的先驗知識，用概率的形式描述目標的運動狀態(tài)[8,10]。本節(jié)在構(gòu)造意圖代價函數(shù)的基礎(chǔ)上，利用貝葉斯理論迭代預測HGRV的機動模式和運動狀態(tài)。

3.1 意圖代價函數(shù)

決定HGRV作戰(zhàn)意圖的主要因素有航向誤差角、彈目距離以及目標重要度[19]。同時，HGRV應(yīng)避免進入一些特定的區(qū)域，如反導系統(tǒng)攔截區(qū)域、易受敵探測或電磁干擾、地緣政治因素不允許通過的區(qū)域。綜合考慮以上因素，定義HGRV意圖代價函數(shù)為：

g(θη|x)exp(-a·σ(x,θη)+b·d(x,θη))·

(5)

式中，θη∈Θ表示打擊目標，φκ∈Φ表示避飛區(qū)；d(·)為彈目距離，σ(·)為航向誤差角；a、b分別表示各因素權(quán)重系數(shù)；Iθ(·)表示打擊目標的重要度；Iφ(x,φκ)表示φκ對代價函數(shù)的影響，定義為：

Iφ(x,φκ)

(6)

式中:c·dmin(x,φκ)表示危險區(qū)半徑，參數(shù)c>1為安全裕度。Iφ(x,φκ)表示當HGRV距離φκ較遠，不考慮避飛區(qū)影響；當HGRV進入危險區(qū)，d(x,φκ)越小，則Iφ(x,φκ)越小，代價越大。下面給出dmin(x,φκ)計算方法。

如圖2所示，令a=1，當HGRV恰好可以以最小轉(zhuǎn)彎半徑r從避飛區(qū)邊緣飛過，此時目標與避飛區(qū)中心O的距離d(x,φκ)定義為dmin(x,φκ)。根據(jù)三角形正弦定理可知：

(7)

式中:HGRV瞬時最小轉(zhuǎn)彎半徑r可根據(jù)狀態(tài)x估算[8]。根據(jù)式(7)可求出dmin(x,φκ)。

3.2 狀態(tài)外推理論

通常認為HGRV會沿代價函數(shù)減小的空域飛行。因此，可根據(jù)代價函數(shù)g(θη|x)定義機動模式j(luò)下狀態(tài)轉(zhuǎn)移概率分布：

p(xi|xi-1,uj)K-1·g(θη|xi)-1

(8)

式中歸一化系數(shù)K為：

(9)

為以一種迭代方式計算目標狀態(tài)的后驗概率密度函數(shù)，計算過程分為模式推斷和狀態(tài)外推。對于模式推斷，首先根據(jù)狀態(tài)分布p(xi+1|x1:i)以及模式的先驗概率p(uj|x1:i-1)推斷模式的后驗概率：

(10)

根據(jù)模式概率Pr(uj|x1:i)對目標狀態(tài)進行一步外推：

(11)

依據(jù)這種思路通過遞推方式更新機動模式和目標狀態(tài)的后驗概率，狀態(tài)的k步外推公式如下：

(12)

3.3 蒙特卡洛算法

由于模式概率p(uj|xi)為非高斯分布，同時，狀態(tài)轉(zhuǎn)移函數(shù)p(xi|xi-1,uj)為非線性的動力學模型，系統(tǒng)通過貝葉斯遞推獲得的狀態(tài)后驗概率函數(shù)p(xi+1|x1:i)無法得到顯式解。作為一種適合于非線性非高斯系統(tǒng)的遞歸貝葉斯濾波方法，蒙特卡洛序貫濾波(粒子濾波)可以利用粒子采樣模擬空間中目標狀態(tài)的后驗概率密度函數(shù)，近似計算目標在空間中概率分布[20]。蒙特卡洛算法的基本步驟如下：

(1) 步驟一 數(shù)據(jù)準備

1) 根據(jù)量測數(shù)據(jù)得到估計值X0，根據(jù)αv和αl預測值構(gòu)造機動模式集U；

2) 參照文獻[21]方法估計HGRV的再入覆蓋范圍，確定可能打擊的目標θη∈Θ；

(2) 步驟二 模式推斷

(3) 步驟三 狀態(tài)外推

3) 根據(jù)式(8)和(9)計算粒子概率。

依據(jù)這種思路對粒子狀態(tài)進行更新和采樣，從而近似得到目標狀態(tài)的概率分布。這種基于粒子概率的描述形式解決了空間中狀態(tài)難以積分運算的問題，提高了算法的適用性。

4 仿真校驗

為驗證方法可行性，設(shè)計以下仿真環(huán)境：1) 飛行器參數(shù)：HGRV模型參考美國洛馬公司CAV-H的基本參數(shù)[22]，采用三自由度動力學方程積分彈道[16]，目標初始速度為4000 m/s，初始高度40 km，攻角隨速度線性變化；2) 量測參數(shù)：定義相控陣雷達采樣間隔為0.1 s，距離量測標準差為500 m，方位角、俯仰角量測標準差均為0.01 rad；3) 算法參數(shù)：預測時間約為340 s，氣動參數(shù)采用ARIMA模型進行時間序列預測，機動模式個數(shù)為5，預測粒子個數(shù)為50；3)仿真場景：HGRV意圖攻擊某目標，采用標準制導法[16]控制傾側(cè)角實現(xiàn)機動轉(zhuǎn)彎并避免飛入特定區(qū)域，飛行時間約為500 s。

實驗1 單目標軌跡預測

圖3～4表示HGRV氣動參數(shù)的估計與預測結(jié)果。不難看出，雖然目標發(fā)生了機動轉(zhuǎn)彎，但αv和αl總體上呈近似線性變化，未發(fā)生劇烈波動，這是因為αv和αl與傾側(cè)角變化無關(guān)。氣動參數(shù)預測RMSE隨時間逐漸增大，但均比預測值小一個量級，在允許范圍內(nèi)?？梢姡鶕?jù)對氣動參數(shù)估計值訓練ARIMA模型對氣動參數(shù)的預測可以較好地逼近真實值，有助于提高軌跡預測準確度。

圖5表示HGRV軌跡的跟蹤與預測粒子散布結(jié)果?？梢钥闯?，為避免進入避飛區(qū)，HGRV不斷改變機動模式飛行。在已知打擊目標的情況下，粒子集根據(jù)飛行意圖進行軌跡預測，當某些粒子距離打擊目標較遠，或者預測粒子進入避飛區(qū)，通過代價函數(shù)計算的粒子權(quán)重逐漸減小，而權(quán)重大的粒子數(shù)量通過重采樣后大幅增加，從而優(yōu)化預測軌跡規(guī)劃結(jié)果。

為對比驗證本文所述方法的有效性，分別采用四種方法對HGRV進行軌跡預測：1) 文獻[4]基于升阻比變化規(guī)律的方法；2) 文獻[5]基于氣動參數(shù)變化規(guī)律的方法；3) 本文所提方法，但不考慮避飛區(qū)；4) 本文所提方法，考慮避飛區(qū)。圖6表示四種方法對比示意圖，其仿真時間分別為：0.24 s、0.27 s、12.83 s、11.75 s。分析上述結(jié)果不難發(fā)現(xiàn)：1) 方法一中升阻比變化規(guī)律未考慮轉(zhuǎn)彎機動，僅對目標縱向狀態(tài)進行外推，預測偏差較大；2) 方法二通過氣動參數(shù)變化規(guī)律辨識目標的機動模式，其軌跡在預測初期與HGRV運動趨勢基本保持一致，當目標改變機動模式后，由于未考慮飛行意圖，預測誤差迅速發(fā)散；3) 方法三和方法四考慮了HGRV的飛行意圖，預測得到的軌跡與目標運動趨勢基本保持一致，雖然預測偏差逐漸增大，但未發(fā)生發(fā)散。同時，由于方法四考慮了避飛區(qū)信息，其預測軌跡更加接近于真實軌跡。4) 從仿真時間來看，本文方法計算復雜度顯著高于傳統(tǒng)算法，這是由于蒙特卡洛方法需要不斷尋優(yōu)最佳機動模式，相對于長達數(shù)百秒的預測軌跡時效而言，仿真時間能夠滿足攔截需求。

實驗2 多目標軌跡預測當HGRV再入覆蓋范圍內(nèi)存在多個可能打擊目標時，根據(jù)方法可以得到多條預測軌跡。為識別HGRV的真實飛行意圖，利用意圖代價函數(shù)計算各目標的概率，進行目標威脅排序。仿真結(jié)果如圖7-8所示，分別表示多目標軌跡預測和意圖推斷示意圖。

從圖7中可以看出，當t=0 s時，根據(jù)算法無法預測意圖為目標1的軌跡，可排除；當t=150 s時，同理可排除目標4；當t=300 s時，對于目標2和3均可以得到預測軌跡，此時應(yīng)同時考慮對兩個目標的威脅。從圖8中可以看出，當時間為0～150 s時，判定目標4威脅最大；隨著HGRV機動轉(zhuǎn)彎，同理依次判定目標1、目標2威脅最大。

以上仿真結(jié)果表明，對于單目標情況，本文方法能夠有效降低目標未來時刻機動不確定對軌跡預測的不利影響。同時，本文方法僅能夠在一定程度上反映目標真實的運動趨勢，無法精確預測軌跡。對于多目標的情況，算法可能存在打擊目標誤判，但隨著跟蹤時間增加，識別準確性會大幅提高。為避免誤判，當多目標概率相差不大時，應(yīng)同時考慮多條預測軌跡，以降低誤判帶來的風險。

5 結(jié) 論

在研究HGRV長期軌跡預測過程中，目標機動是難以回避的問題。本文在傳統(tǒng)狀態(tài)外推預測思路的基礎(chǔ)上，提出基于意圖推斷的貝葉斯軌跡預測方法。方法根據(jù)貝葉斯理論迭代求解了HGRV的運動狀態(tài)和機動模式，利用蒙特卡洛采樣近似計算了HGRV的概率分布，避免了復雜空域概率密度積分困難。

仿真對比表明，由于升力參數(shù)和阻力參數(shù)變化較小，利用ARIMA模型對其進行參數(shù)預測取得了預期效果。相對傳統(tǒng)方法，結(jié)合意圖推斷的軌跡預測方法有助于減小目標未知機動對軌跡預測的不利影響。對于多目標的情況，應(yīng)同時考慮其打擊的可能性，降低誤判帶來的風險。