基于意圖推斷的高超聲速滑翔目標(biāo)貝葉斯軌跡預(yù)測

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

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

0 引 言

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

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

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

1 問題描述

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

2 預(yù)測模型

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

2.1 動力學(xué)模型

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

(1)

(2)

式中:r為目標(biāo)瞬時地心距；B為雷達(dá)站地理緯度；Re為地球半徑。

2.2 機(jī)動模式集

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

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

(3)

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

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

(4)

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

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

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

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

3 貝葉斯軌跡預(yù)測

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

3.1 意圖代價函數(shù)

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

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

(5)

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

Iφ(x,φκ)

(6)

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

如圖2所示，令a=1，當(dāng)HGRV恰好可以以最小轉(zhuǎn)彎半徑r從避飛區(qū)邊緣飛過，此時目標(biāo)與避飛區(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)外推理論

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

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

(8)

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

(9)

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

(10)

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

(11)

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

(12)

3.3 蒙特卡洛算法

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

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

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

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

(2) 步驟二 模式推斷

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

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

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

4 仿真校驗

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

實驗1 單目標(biāo)軌跡預(yù)測

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

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

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

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

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

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

5 結(jié) 論

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

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