臨近空間升力式再入飛行器跳躍空域預(yù)測(cè)

喻晨龍， 譚賢四， 曲智國， 謝 非， 李陸軍

(1. 空軍預(yù)警學(xué)院研究生管理大隊(duì)，湖北 武漢 430019； 2. 空軍預(yù)警學(xué)院防空預(yù)警裝備系， 湖北 武漢430019；3. 長江勘測(cè)規(guī)劃設(shè)計(jì)研究院， 湖北 武漢 430019)

臨近空間升力式再入飛行器的飛行高度高、飛行速度快、突防能力強(qiáng)、打擊距離遠(yuǎn)，具有升力式氣動(dòng)布局，其由火箭助推至高空后，依靠地球重力再入大氣層，可長距離、無動(dòng)力滑翔飛行，典型的有HTV-2和AHW兩種[1-3]。這類再入式飛行器在氣動(dòng)外形、運(yùn)行環(huán)境、飛行受力、控制方式和機(jī)動(dòng)模式等方面與彈道導(dǎo)彈存在明顯差異，給預(yù)警探測(cè)系統(tǒng)帶來了極大的挑戰(zhàn)[4-8]。針對(duì)這一殺手锏武器，目前可行的探測(cè)手段是以反導(dǎo)預(yù)警系統(tǒng)為依托，以米波雷達(dá)、凝視雷達(dá)等新型裝備為補(bǔ)充進(jìn)行協(xié)同探測(cè)，從而盡早盡快發(fā)現(xiàn)目標(biāo)并維持較高的跟蹤精度，為后端的攔截系統(tǒng)提供連續(xù)穩(wěn)定的信息支撐。

通常，遠(yuǎn)程預(yù)警相控陣?yán)走_(dá)靠近邊境線前沿部署，探測(cè)距離可達(dá)5 000 km，然而受地球曲率的限制，在探測(cè)臨近空間目標(biāo)時(shí)探測(cè)距離最遠(yuǎn)只有1 300 km。臨近空間升力式再入飛行器在縱向上存在跳躍機(jī)動(dòng)，在橫向上存在側(cè)向偏移，其中：只有當(dāng)目標(biāo)處于雷達(dá)探測(cè)范圍內(nèi)時(shí)，才能保持連續(xù)穩(wěn)定跟蹤；當(dāng)目標(biāo)飛出探測(cè)范圍，即向下機(jī)動(dòng)至地平線以下時(shí)會(huì)使情報(bào)中斷。若能利用目標(biāo)在雷達(dá)水平面波束層中飛行時(shí)獲得的有限觀測(cè)信息，對(duì)目標(biāo)的運(yùn)動(dòng)包絡(luò)和可能的跳躍空域進(jìn)行預(yù)測(cè)，則可輔助雷達(dá)在遠(yuǎn)程預(yù)警相控陣?yán)走_(dá)丟失目標(biāo)后補(bǔ)充照射，從而保證情報(bào)的連續(xù)。

升力式再入飛行器的彈道包括平衡滑翔和跳躍滑翔[9]。研究者們針對(duì)目標(biāo)的運(yùn)動(dòng)參數(shù)、攻角和傾側(cè)角的變化規(guī)律進(jìn)行了深入分析，如：李廣華等[10]認(rèn)為，在平衡滑翔狀態(tài)下目標(biāo)的速度與高度、速度與速度傾角之間存在解析關(guān)系，且隨著傾側(cè)角的變化而變化；文獻(xiàn)[11-12]的作者通過分析氣動(dòng)力加速度對(duì)目標(biāo)質(zhì)心運(yùn)動(dòng)和繞質(zhì)心運(yùn)動(dòng)的影響，說明氣動(dòng)力改變是目標(biāo)跳躍機(jī)動(dòng)的主要原因；李廣華[13]認(rèn)為目標(biāo)的運(yùn)動(dòng)狀態(tài)參數(shù)反映著目標(biāo)的氣動(dòng)力，只有當(dāng)參數(shù)滿足特定數(shù)學(xué)關(guān)系時(shí)，才能維持平衡滑翔軌跡，當(dāng)偏離平衡點(diǎn)時(shí)就會(huì)產(chǎn)生跳躍軌跡；王路等[14]通過目標(biāo)的三自由度彈道方程計(jì)算了目標(biāo)的攻角和最大升阻比。然而，上述研究未利用運(yùn)動(dòng)參數(shù)、攻角和傾側(cè)角的變化規(guī)律進(jìn)行彈道預(yù)測(cè)。鑒于此，筆者在分析目標(biāo)的運(yùn)動(dòng)特性的基礎(chǔ)上，估算目標(biāo)的運(yùn)動(dòng)參數(shù)、攻角和傾側(cè)角，合理假設(shè)攻角和傾側(cè)角的變化規(guī)律，最后通過數(shù)值仿真方法預(yù)測(cè)目標(biāo)飛離雷達(dá)探測(cè)范圍后可能的跳躍空域。

1 模型建立

1.1 飛行器的三自由度運(yùn)動(dòng)模型

為了描述目標(biāo)的運(yùn)動(dòng)過程和雷達(dá)的探測(cè)過程，定義地心(Earth Canter，EC)坐標(biāo)系、飛行器東北天(Vehicle East-North-Up，VENU)坐標(biāo)系和雷達(dá)站東北天(Radar East-North-Up, RENU)坐標(biāo)系3個(gè)相關(guān)的參考坐標(biāo)系[14]。這3個(gè)坐標(biāo)系各自包括1個(gè)直角坐標(biāo)系和1個(gè)球坐標(biāo)系。

假設(shè)地球?yàn)榫鶆驁A球，不考慮自轉(zhuǎn)影響，則飛行器在無動(dòng)力滑翔過程中滿足EC球坐標(biāo)系和VENU直角坐標(biāo)系的三自由度運(yùn)動(dòng)方程[8，14],即

(1)

(2)

式中：r(t)、χ(t)、φ(t)分別為飛行器在EC球坐標(biāo)系下的地心距、經(jīng)度和緯度；v(t)、γ(t)、ψ(t)分別為飛行器在VENU直角坐標(biāo)系下的速度、航跡傾角和航向角；m為飛行器的質(zhì)量；g為重力加速度；σ(t)為傾側(cè)角；D(t)、L(t)分別為飛行過程中所受阻力和升力,

(3)

(4)

其中CD(t)為阻力系數(shù)，CL(t)為升力系數(shù)，ρ(t)為大氣密度，S為飛行器氣動(dòng)參考面積。

文獻(xiàn)[15-18]作者認(rèn)為：當(dāng)飛行器的速度>8 Ma時(shí)，升力系數(shù)可簡化為攻角α(t)的線性函數(shù)，阻力系數(shù)可簡化為α(t)的二次函數(shù)。對(duì)氣動(dòng)參數(shù)插值擬合后可得

CD(t)=0.044+0.08α(t)+2.2α2(t),

(5)

CL(t)=-0.151 7+3.151α(t)。

(6)

大氣密度滿足指數(shù)模型，即

ρ(t)=ρ0exp[-(r(t)-Re)/hs]，

(7)

式中：ρ0=1.752 kg/m3，為參考大氣密度；hs=6.7 km，為參考高度；Re=6 371 km，為地球半徑。

1.2 雷達(dá)的觀測(cè)模型

飛行器穿越水平波束屏的過程中，雷達(dá)每間隔一段時(shí)間Ts(即搜索間隔時(shí)間)照射一次目標(biāo)，得到了N組觀測(cè)值{z1,z2,…,zN}，這些觀測(cè)值是在RENU球坐標(biāo)系下獲取的，記為

zi=[?r(ti)θr(ti)φr(ti)]T，

(8)

式中：?r(ti)、θr(ti)、φr(ti) 分別為第i次照射時(shí)的目標(biāo)徑向距離、方位角、俯仰角，其中ti為第i次照射時(shí)對(duì)應(yīng)的時(shí)刻。將觀測(cè)值在RENU直角坐標(biāo)系中表示，有

(9)

通過式(9)位置差分可得到速度信息Vr(ti)=[vrx(ti)vry(ti)vrz(ti)]T。

2 參數(shù)估計(jì)

防御方預(yù)警系統(tǒng)可直接獲取的信息僅包括目標(biāo)在穿越水平波束屏過程中的觀測(cè)信息，即RENU直角坐標(biāo)系下的位置和速度。而目標(biāo)的面質(zhì)比、控制量，以及在EC球坐標(biāo)系下的地理位置和VENU直角坐標(biāo)系下的速度是未知的。

由式(1)-(7)可知：目標(biāo)的彈道軌跡取決于運(yùn)動(dòng)狀態(tài)初值和攻角/傾側(cè)角的變化規(guī)律。為了實(shí)現(xiàn)軌跡預(yù)測(cè)，首先把觀測(cè)信息由RENU直角坐標(biāo)系轉(zhuǎn)換到EC球坐標(biāo)系和VENU直角坐標(biāo)系，得到目標(biāo)飛離水平波束屏的運(yùn)動(dòng)狀態(tài)，然后結(jié)合式(1)-(7)估算目標(biāo)的攻角/傾側(cè)角和面質(zhì)比。

2.1 運(yùn)動(dòng)參數(shù)初值估計(jì)

1) 將位置觀測(cè)信息從RENU直角坐標(biāo)系轉(zhuǎn)換到EC球坐標(biāo)系。首先從RENU直角坐標(biāo)系轉(zhuǎn)換到EC直角坐標(biāo)系，然后從EC直角坐標(biāo)系轉(zhuǎn)換到EC球坐標(biāo)系。

假設(shè)雷達(dá)站OR所在地的經(jīng)度為J，緯度為W，海拔高度為H。要從RENU直角坐標(biāo)系變換到EC直角坐標(biāo)系，根據(jù)歐拉角關(guān)系，首先把RENU直角坐標(biāo)系原點(diǎn)先沿地心矢移動(dòng)Re+H到地心O，然后繞x軸順時(shí)針旋轉(zhuǎn)π/2-W，最后繞z軸順時(shí)針旋轉(zhuǎn)π/2+J。即

(10)

式中：Xr(ti)=[xr(ti)yr(ti)zr(ti)]T，為RNEU直角坐標(biāo)系下的位置；Xe(ti)=[xe(ti)ye(ti)ze(ti)]T，為EC直角坐標(biāo)系下的位置；

(11)

為從RENU直角坐標(biāo)系到EC直角坐標(biāo)系的坐標(biāo)變換矩陣，其中Tx、Tz分別為繞x、z軸的基元旋轉(zhuǎn)矩陣。

將位置觀測(cè)信息從EC直角坐標(biāo)系轉(zhuǎn)換到EC球坐標(biāo)系，即

(12)

假設(shè)在第k次照射后飛行器飛離水平波束屏，雷達(dá)測(cè)量得到目標(biāo)的質(zhì)心位置為[?r(tk)θr(tk)φr(tk)]T(其中tk為第k次照射時(shí)對(duì)應(yīng)的時(shí)刻)，這是一個(gè)RENU球坐標(biāo)，通過坐標(biāo)變換得到此時(shí)質(zhì)心在EC球坐標(biāo)系下位置為[re(tk)χe(tk)φe(tk)]T。

2) 將速度信息從RENU直角坐標(biāo)系轉(zhuǎn)換到VENU直角坐標(biāo)系。先考慮位置坐標(biāo)的變換，將位置坐標(biāo)從RENU直角坐標(biāo)系轉(zhuǎn)換到EC直角坐標(biāo)系，再從EC直角坐標(biāo)系轉(zhuǎn)換到VENU直角坐標(biāo)系。由于從EC直角坐標(biāo)系轉(zhuǎn)換到VENU直角坐標(biāo)系，與從RENU直角坐標(biāo)系轉(zhuǎn)換到EC直角坐標(biāo)系的變換關(guān)系剛好相反，因此可類比式(10)中的方法轉(zhuǎn)換。特別需注意的是，此時(shí)坐標(biāo)系的原點(diǎn)變?yōu)榱孙w行器的質(zhì)心。

要從EC直角坐標(biāo)系轉(zhuǎn)換到VENU直角坐標(biāo)系，根據(jù)歐拉角關(guān)系，首先把EC直角坐標(biāo)系先繞z軸逆時(shí)針旋轉(zhuǎn)π/2+χe(tk)，然后繞x軸逆時(shí)針旋轉(zhuǎn)π/2-φe(tk)，最后坐標(biāo)原點(diǎn)由地心O沿地心矢移動(dòng)re(tk)到飛行器的質(zhì)心。即

(13)

式中：Xe(ti)=[xe(ti)ye(ti)ze(ti)]T和Xv(ti)=[xv(ti)yv(ti)zv(ti)]T分別表示第i次照射時(shí)目標(biāo)在EC直角坐標(biāo)系和VENU直角坐標(biāo)系下的位置；

(14)

為從EC直角坐標(biāo)系到VENU直角坐標(biāo)系的坐標(biāo)變換矩陣。

那么，從RENU直角坐標(biāo)系到VENU直角坐標(biāo)系的位置坐標(biāo)變換關(guān)系為

(15)

因此，從RENU直角坐標(biāo)系下的速度矢量[vrx(ti)vry(ti)vrz(ti)]T轉(zhuǎn)換到VENU直角坐標(biāo)系下的速度矢量[vvx(ti)vvy(ti)vvz(ti)]T，其變換關(guān)系為

(16)

由于VNEU直角坐標(biāo)系下的速度、航跡傾角與航向角有如下關(guān)系:

(17)

因此，根據(jù)RENU直角坐標(biāo)系下的一組量測(cè)信息差分得到飛行器飛離水平波束屏的速度[vrx(tk)vry(tk)vrz(tk)]T，然后根據(jù)坐標(biāo)轉(zhuǎn)換關(guān)系把各軸上的分速度變換到VENU直角坐標(biāo)系下，最后得到飛行器的速度vv(tk)、航跡傾角γv(tk)和航向角ψv(tk)。

2.2 攻角和傾側(cè)角估計(jì)

式(1)、(2)描述了目標(biāo)的地理位置和速度的迭代關(guān)系，式(1)-(7)隱含了目標(biāo)控制量的解算方法。因此，將式(2)改寫為[20]

(18)

(19)

(20)

式中：ρ(ti)=ρ0exp[-(re(ti)-Re)/hs]。

(21)

(22)

(23)

式中：K(ti)為目標(biāo)的升阻比。

將式(21)代入式(5)、(6)，可以解算出目標(biāo)第i次照射時(shí)的攻角α(ti)、升力系數(shù)CL(ti)和阻力系數(shù)CD(ti)。

3 跳躍空域預(yù)判

當(dāng)初值確定后，飛行器的彈道由攻角和傾側(cè)角決定，攻角和傾側(cè)角并不是隨機(jī)變化的。其中：1) 攻角直接決定飛行過程中的氣動(dòng)特性，再入過程中若調(diào)整頻繁且幅值較大，將嚴(yán)重影響飛行器的熱防護(hù)系統(tǒng)，因此攻角的變化要盡量簡單，一般選擇分段線性函數(shù)，大小通?？刂茷?°～ 20°，此外攻角與升阻比緊密相關(guān)，最常用的制導(dǎo)律是以最大升阻比飛行;2) 傾側(cè)角是主控制量，其確定升力在縱向和側(cè)向上的分量，并決定飛行器橫向機(jī)動(dòng)程度，改變傾側(cè)角是為了規(guī)避禁飛區(qū)和突防，傾側(cè)角最常用的制導(dǎo)律是方位角誤差走廊，其幅值不變、符號(hào)翻轉(zhuǎn)，幅值通常控制在±40°以內(nèi)。

因此，升力式再入飛行器在遠(yuǎn)端運(yùn)動(dòng)時(shí)，為了維持最小能耗和最大升阻比飛行，并且實(shí)現(xiàn)有效突防，對(duì)攻角和傾側(cè)角的變化規(guī)律可做如下假設(shè)：

(24)

(25)

式中：α(tk)為飛行器飛出水平玻束屏?xí)r估計(jì)的攻角；αmaxK(tk)為最大升阻比對(duì)應(yīng)的攻角；σ(tk)為飛行器飛出水平玻束屏?xí)r估計(jì)的傾側(cè)角；t為飛行時(shí)間；T為攻角改變的時(shí)間；T1為傾側(cè)角改變的時(shí)間。

將飛行器飛離水平波束屏?xí)r的運(yùn)動(dòng)參數(shù)、攻角和傾側(cè)角的預(yù)測(cè)值代入式(1)-(7)進(jìn)行數(shù)值仿真，隨著迭代的進(jìn)行，當(dāng)式(2)中第1次出現(xiàn)γ(t)=0時(shí)停止。

4 仿真驗(yàn)證

假設(shè)某型遠(yuǎn)程預(yù)警相控陣?yán)走_(dá)部署地點(diǎn)在東經(jīng)110°、北緯10°、海拔高度0 km，波束寬度1.2°，沿著正東設(shè)置單層水平波束屏，量測(cè)噪聲服從高斯分布，測(cè)量誤差??=100 m，?θ=?φ=0.001 5°，搜索間隔時(shí)間Ts=1 s。飛行器穿越水平波束屏的過程中雷達(dá)在RENU球坐標(biāo)系內(nèi)得到了14組觀測(cè)數(shù)據(jù)，如表1所示。

表1 RENU球坐標(biāo)系內(nèi)觀測(cè)數(shù)據(jù)

按照2.1節(jié)的估算方法對(duì)觀測(cè)數(shù)據(jù)進(jìn)行坐標(biāo)轉(zhuǎn)換，計(jì)算得到目標(biāo)飛離水平波束屏?xí)r的地理位置和速度分別為：re(tk)=6 425 km，χe(tk)=117.472 8°，φe(tk)=9.889 5°；vv(tk)=19.850 9 Ma，γv(tk)=-1.868 6°，ψv(tk)=92.029 4°。進(jìn)一步，根據(jù)2.2節(jié)的估算方法，計(jì)算得到目標(biāo)的面質(zhì)比為S/m=0.002 5，在飛離水平波束屏?xí)rα(tk)≈5°，σ(tk)≈40°。

飛行器飛離水平波束屏后，其攻角和傾側(cè)角按照式(24)、(25)的規(guī)律變化，即攻角和傾側(cè)角先是維持不變，再分別在時(shí)刻T和T1發(fā)生改變，攻角由飛離時(shí)的攻角α(tk) 變化為最大升阻比攻角αmaxK(tk) ，傾側(cè)角幅值不變、符號(hào)翻轉(zhuǎn)。這時(shí)，隨著迭代的進(jìn)行，軌跡發(fā)生較大變化。由于攻角和傾側(cè)角改變的時(shí)刻是無法預(yù)知的，當(dāng)式(2)中第1次出現(xiàn)γ(t)=0時(shí)迭代停止，此時(shí)的軌跡曲線就圍成了一個(gè)空域范圍Ω，在EC球坐標(biāo)系下的目標(biāo)跳躍空域范圍如圖1所示。

5 結(jié)論

臨近空間武器的興起將打破常規(guī)作戰(zhàn)規(guī)則，在探測(cè)這類目標(biāo)時(shí)，遠(yuǎn)程預(yù)警雷達(dá)因受到地球地平線遮擋而丟失目標(biāo)，本文基于目標(biāo)彈道特性，利用有限的觀測(cè)數(shù)據(jù)對(duì)目標(biāo)可能的跳躍空域進(jìn)行預(yù)判，能提前指示預(yù)警系統(tǒng)調(diào)度雷達(dá)資源探測(cè)，具有較強(qiáng)的軍事應(yīng)用價(jià)值。

參考文獻(xiàn)：

[1] PAUL L M,VINCENT L R,LUAT T N,et al.NASA hypersonic fight demonstration overview,status,and future plans[J].Acta astronautica,2004,55(7):619-630.

[2] 李益翔.美國高超聲速飛行器發(fā)展歷程研究[D].哈爾濱：哈爾濱工業(yè)大學(xué),2016：18-32.

[3] STEVEN H W,JEFFREY S,DALE S,et al.The DAPPA/AF falcon program:the hypersonic technology vehicle #2(HTV-2) flight

demonstration phase[C]∥15th AIAA International Space Planes and Hypersonic Systems and Technologies Conference.28 April-1 2008,Dayton,Ohio:2008-2539.

[4] RANDALL T V,LAERENCE D H,CHARLES R,M,et al.X-43A Hypersonic vehicle technology development[J].Acta Astronautica,2006(59):181-191.

[5] STEVEN W,MING T,SUE M,et al.Falcon HTV-3X-A reusable hypersonic test bed[C]∥15th AIAA International Space Planes and Hypersonic Systems and Technologies Conference.28 April-1 May,2008,Dayton,Ohio：2008-2544.

[6] 秦雷，李君龍，周荻.國外臨近空間高速滑翔彈頭運(yùn)動(dòng)特征研究[J].戰(zhàn)術(shù)導(dǎo)彈技術(shù),2015(1):16-20.

[7] 李君龍，秦雷，謝曉瑛.臨近空間非彈道式機(jī)動(dòng)模式與跟蹤濾波問題[J].現(xiàn)代防御技術(shù),2017,45(1):99-106.

[8] 李慧峰.高超聲速飛行器制導(dǎo)與控制技術(shù)[M].北京：中國宇航出版社,2012:10.

[9] 陳小慶，候中喜，劉建霞.高超聲速滑翔飛行器彈道特性分析[J].導(dǎo)彈與航天運(yùn)載技術(shù),2011(2):5-9.

[10] 李廣華，張洪波，湯建國.高超聲速滑翔飛行器典型彈道特性分析[J].宇航學(xué)報(bào)，2015,36(4):397-403.

[11] OLIVIER D M,ROBERT H B.Tracking and identification of a maneuvering reentry vehicle[C]∥AIAA Guidance,Navigation,and Control Conference and Exhibit.11-14 August 2003,Austin,Texas,USA,2003:2003-5446.

[12] 張凱，熊家軍，韓春耀，等.一種基于氣動(dòng)力模型的高超聲速滑翔目標(biāo)跟蹤算法[J].宇航學(xué)報(bào),2017,38(2):123-130.

[13] 李廣華.高超聲速滑翔飛行器運(yùn)動(dòng)特性分析及彈道跟蹤預(yù)報(bào)方法研究[D].長沙:國防科學(xué)技術(shù)大學(xué),2016：17-49.

[14] 王路，邢清華，毛藝帆.基于升阻比變化規(guī)律的再入高超聲速滑翔飛行器軌跡預(yù)測(cè)算法[J].系統(tǒng)工程與電子技術(shù),2015,37(10):2335-2339.

[15] PHILIPS T H.A Common Aero Vehicle(CAV) model,description,and employment guide[R].USA:Schafer Corporation,2003.

[16] 趙漢元.飛行器再入動(dòng)力學(xué)與制導(dǎo)[M].長沙:國防科技大學(xué)出版社,1997：48-97.

[17] 黃長強(qiáng)，國海峰，丁達(dá)理.高超聲速滑翔飛行器軌跡優(yōu)化與制導(dǎo)綜述[J].宇航學(xué)報(bào),2014,35(4):369-379.

[18] 王璐璐，秦玉亮，王宏強(qiáng)，等.飛行器升阻比估計(jì)誤差分析方法研究[J].雷達(dá)科學(xué)與技術(shù),2012,10(2):175-179.