吳 楠，王 鋒，趙 敏，孟凡坤

(1.中國(guó)人民解放軍戰(zhàn)略支援部隊(duì)信息工程大學(xué) 數(shù)據(jù)與目標(biāo)工程學(xué)院，河南 鄭州 450001;2.中國(guó)人民解放軍95894部隊(duì)，北京 102211)

高超聲速滑翔再入飛行器[1-2](Hypersonic Glide Reentry Vehicle,HGRV)的可達(dá)區(qū)是指在某些約束條件下，飛行器在地球表面上可達(dá)范圍的邊界曲線[3]?？蛇_(dá)區(qū)作為飛行器縱橫向機(jī)動(dòng)能力和打擊區(qū)域的反映，其計(jì)算具有重要意義。

可達(dá)區(qū)計(jì)算的內(nèi)容包括兩個(gè)方面，一是最大橫程的計(jì)算方法，二是邊界的獲取方式。其中前者的關(guān)鍵是獲得傾側(cè)角的變化規(guī)律，目前有三種方法：①通過(guò)偽譜[4-5]或粒子群[6]等數(shù)值優(yōu)化方法獲得傾側(cè)角的變化規(guī)律；②利用最優(yōu)化原理[7]或再入走廊邊界[8]推出傾側(cè)角的控制律表達(dá)式；③將傾側(cè)角設(shè)置為常值[9]。這三種方法的精度、運(yùn)算量和依賴的先驗(yàn)信息量是逐漸減小的。可達(dá)區(qū)邊界的獲取目前有兩種方法：①遍歷法，即通過(guò)計(jì)算不同縱程或傾側(cè)角條件下的最大橫程彈道，將末段點(diǎn)連接構(gòu)成邊界；②橢圓近似法，即將可達(dá)區(qū)近似為橢圓形，利用最大縱程和最大橫程三個(gè)末端點(diǎn)計(jì)算橢圓參數(shù)形成邊界[10]。相比之下遍歷法計(jì)算的邊界與實(shí)際更為符合，但帶來(lái)的問(wèn)題是運(yùn)算量急劇增加。

傳統(tǒng)的可達(dá)區(qū)計(jì)算方法通常為上述方法的不同組合，在計(jì)算精度和運(yùn)算量上具有較大差異，但均是從飛行器設(shè)計(jì)方的角度，即在已知詳細(xì)的參數(shù)和約束條件下進(jìn)行構(gòu)建，重視準(zhǔn)確性而忽略時(shí)效性，通常用于飛行器的彈道設(shè)計(jì)與性能分析。然而對(duì)于防御方，對(duì)來(lái)襲的HGRV進(jìn)行預(yù)警和威脅評(píng)估也需要對(duì)其可達(dá)區(qū)進(jìn)行預(yù)測(cè)，而傳統(tǒng)的可達(dá)區(qū)計(jì)算方法不太適用，原因?yàn)椋阂环矫妫藭r(shí)HGRV作為非合作目標(biāo)，相關(guān)參數(shù)和約束是未知的；另一方面，可達(dá)區(qū)計(jì)算作為針對(duì)HGRV的“預(yù)警—探測(cè)—預(yù)判—攔截”中重要一環(huán)，要求需具有準(zhǔn)實(shí)時(shí)性，傳統(tǒng)方法具有依賴先驗(yàn)信息偏多、耗時(shí)較長(zhǎng)的缺陷。

本文在僅已知目標(biāo)當(dāng)前時(shí)刻的位置、速度和最大升阻比參數(shù)(可基于雷達(dá)探測(cè)數(shù)據(jù)通過(guò)實(shí)時(shí)彈道估計(jì)獲得)條件下，對(duì)文獻(xiàn)[11]中的傾側(cè)角最優(yōu)控制律進(jìn)行簡(jiǎn)化和改進(jìn)，并引入平衡滑翔和最大升阻比滑翔假設(shè)，建立簡(jiǎn)化的飛行程序控制模型，分別通過(guò)一次數(shù)值積分獲得最大縱程和橫程彈道，基于橢圓近似法利用三個(gè)末端點(diǎn)構(gòu)建可達(dá)區(qū)橢圓邊界，以實(shí)現(xiàn)對(duì)可達(dá)區(qū)的準(zhǔn)實(shí)時(shí)計(jì)算。

1 彈道計(jì)算方程的簡(jiǎn)化

HGRV的三自由度彈道計(jì)算方程通常如式(1)所示。

(1)

當(dāng)飛行器飛至中末段或飛行器升阻比有限，而使得剩余飛行時(shí)間較短時(shí)，可以忽略地球自轉(zhuǎn)，式(1)可簡(jiǎn)化為

(2)

2 基于最大升阻比平衡滑翔的最大縱程

飛行器最大縱程的計(jì)算可以利用三個(gè)條件：①傾側(cè)角ν=0，即飛行器在縱向平面運(yùn)動(dòng)，不進(jìn)行橫向滾轉(zhuǎn)機(jī)動(dòng)；②平衡滑翔(Equilibrium Glide,EG)假設(shè)，彈道平穩(wěn)便于控制，且過(guò)程約束容易滿足；③最大升阻比滑翔的近似最優(yōu)性，與最大縱程優(yōu)化結(jié)果誤差小于2%。

平衡滑翔[12]是指飛行器升力與重力達(dá)到某種平衡，從而實(shí)現(xiàn)平穩(wěn)滑翔的狀態(tài)，其優(yōu)點(diǎn)是利用該假設(shè)可顯著降低彈道計(jì)算的復(fù)雜度，提高運(yùn)算速度，且計(jì)算所得的彈道通常滿足飛行過(guò)程中的動(dòng)壓、熱流和過(guò)載等約束條件，經(jīng)常用于再入飛行器滑翔彈道的設(shè)計(jì)與分析。

平衡滑翔的條件為

(3)

對(duì)于最大縱程，ν=0，并綜合式(2)和式(3)可得平衡滑翔條件下的升力表達(dá)式為

(4)

可以看出，平衡滑翔條件下的升力隨著目標(biāo)的速度、高度和速度傾角變化而變化，但不一定與最大升阻比對(duì)應(yīng)的升力相等，說(shuō)明平衡滑翔條件可使彈道較為平穩(wěn)，但射程不一定最大。為獲得最大縱程，本文提出一種“最大升阻比平衡滑翔”的虛擬狀態(tài)，該狀態(tài)的升力由平衡滑翔條件計(jì)算，而氣動(dòng)阻力則由升力和最大升阻比常數(shù)表示

(5)

由于Kmax為目標(biāo)升阻比的極大值，因此Dmin為實(shí)際氣動(dòng)阻力的極小值，“最大升阻比平衡滑翔”的本質(zhì)就是恒定氣動(dòng)阻力極小值條件下的平衡滑翔，從而融合平衡滑翔和最大升阻比滑翔的優(yōu)勢(shì)，達(dá)到既滿足過(guò)程約束又近似最大縱程的目的。

將式(5)代入式(2)可得

(6)

將式(4)和式(6)代入式(2)，可得基于最大升阻比平衡滑翔條件下的最大縱程簡(jiǎn)化計(jì)算的積分方程

(7)

(8)

式中，S和q分別為目標(biāo)的截面積和動(dòng)壓。

給定某一終端速度約束Vf，便可以利用式(7)通過(guò)數(shù)值積分算法計(jì)算獲得最大縱程的終點(diǎn)坐標(biāo)(λ1,φ1)。

3 基于最優(yōu)解擬合公式的最大橫程計(jì)算

基于最優(yōu)化原理，文獻(xiàn)[11]推導(dǎo)出橫程的埃格斯公式

(9)

可以看出，橫程僅為最大升阻比Kmax和傾側(cè)角ν的函數(shù)，如果忽略掉級(jí)數(shù)項(xiàng)，使橫程取得極大值的傾側(cè)角為ν=45°，這也是工程中常用ν=45°的常值傾側(cè)角控制來(lái)獲得近似最大橫程的原因。當(dāng)最大升阻比較大時(shí)，忽略級(jí)數(shù)項(xiàng)會(huì)帶來(lái)較大誤差，實(shí)際上最優(yōu)的傾側(cè)角控制律中，傾側(cè)角是時(shí)變的，且與最大升阻比有關(guān)。此時(shí)最優(yōu)傾側(cè)角的計(jì)算公式可表示為

(10)

式中：βC為橫程角，初始為0，隨著橫程增加而增大；σ-σ0表征速度矢量與飛行縱向平面的夾角，初始為0，逐漸增加至約90°。整體來(lái)看，傾側(cè)角隨著速度矢量與縱向平面夾角的增加而減小，且在速度矢量與縱向平面垂直前減小至0，防止航程回旋而導(dǎo)致橫程減小。而最大升阻比則決定了初始傾側(cè)角的大小，最大升阻比越大，初始傾側(cè)角也越大，說(shuō)明飛行器用于橫向機(jī)動(dòng)的偏轉(zhuǎn)能力越大。

當(dāng)傾側(cè)角不為0時(shí)，將平衡滑翔約束式(3)代入式(2)可得此時(shí)的升力表達(dá)式

(11)

(12)

仍然基于“最大升阻比平衡滑翔”，可得傾側(cè)角不為0時(shí)氣動(dòng)阻力用最大升阻比表示的表達(dá)式

(13)

將式(11)和式(13)代入式(2)可得最大橫程計(jì)算的積分方程

(14)

式(10)和式(14)構(gòu)成了閉合的最大橫程彈道計(jì)算方程，方程中仍只有一個(gè)未知參數(shù)即最大升阻比，采用數(shù)值積分算法，即可分別計(jì)算左向和右向兩條最大橫程彈道，彈道終點(diǎn)坐標(biāo)分別為(λ2,φ2)和(λ3,φ3)。

4 可達(dá)區(qū)橢圓計(jì)算

根據(jù)球面幾何，已知球面上三點(diǎn)計(jì)算球面上的橢圓區(qū)域是非常復(fù)雜的，考慮到本文研究的飛行器可達(dá)區(qū)域遠(yuǎn)小于軌道飛行器的可達(dá)區(qū)，因此假設(shè)飛行器的可達(dá)區(qū)近似在以落點(diǎn)為中心、經(jīng)度軸為橫軸、緯度軸為縱軸的二維平面區(qū)域。

在該經(jīng)緯度二維平面內(nèi)，定義可達(dá)區(qū)為以(λ2,φ2)和(λ3,φ3)間的線段為短軸、(λ1,φ1)到(λ2,φ2)和(λ3,φ3)連線中點(diǎn)的距離為半長(zhǎng)軸的半橢圓區(qū)域。

橢圓中心的坐標(biāo)為

(15)

橢圓的長(zhǎng)半軸和短半軸為

(16)

該橢圓相當(dāng)于坐標(biāo)在原點(diǎn)的標(biāo)準(zhǔn)橢圓先旋轉(zhuǎn)Φ角度，然后再將中心平移至(λ0,φ0)。Φ的計(jì)算公式為

(17)

可達(dá)區(qū)橢圓滿足公式

(18)

5 仿真分析

假設(shè)地面雷達(dá)對(duì)某來(lái)襲HGRV進(jìn)行一段穩(wěn)定跟蹤后，通過(guò)濾波和實(shí)時(shí)彈道估計(jì)，獲得其當(dāng)前狀態(tài)參數(shù)估計(jì)值如表1所示。

表1 目標(biāo)當(dāng)前時(shí)刻狀態(tài)參數(shù)估計(jì)值Tab.1 State estimation of target at current time

首先將當(dāng)前狀態(tài)參數(shù)估計(jì)轉(zhuǎn)化為本文算法所需的初值，即

[r,φ,λ,V,Θ,σ]=f([X,Y,Z,VX,VY,VZ])

(19)

利用式(7)和式(14)進(jìn)行數(shù)值積分，根據(jù)目標(biāo)打擊通常對(duì)終端速度的要求，定Vf=1 800 m/s為積分終止條件，便可獲得縱程終點(diǎn)坐標(biāo)和兩個(gè)橫程終點(diǎn)坐標(biāo)，進(jìn)而計(jì)算可達(dá)區(qū)橢圓，本文算法計(jì)算的HGRV最大縱程、橫程以及可達(dá)區(qū)橢圓結(jié)果如表2和圖1所示。

表2 可達(dá)區(qū)橢圓參數(shù)Tab.2 Parameters of footprint ellipse 單位：(°)

圖1 最大縱程和橫程地面航跡與可達(dá)區(qū)橢圓邊界Fig.1 Ground track of maximum downrange and crossrange and elliptical boundary of footprint

從結(jié)果可以看出，對(duì)于本文算例中的HGRV(初始速度約為3 800 m/s，最大升阻比為3)，其最大縱程約1 960 km，最大橫程約為650 km，可達(dá)區(qū)面積(橢圓的前半段)仍可達(dá)9×105km2。

為檢驗(yàn)本文算法的準(zhǔn)確性和有效性，分別采用工程中常用的數(shù)值優(yōu)化方法和常傾側(cè)角方法[14]計(jì)算可達(dá)區(qū)進(jìn)行比較。其中，數(shù)值優(yōu)化方法計(jì)算精度高(但需要利用目標(biāo)充分的總體參數(shù)信息)，可作為結(jié)果準(zhǔn)確度驗(yàn)證的依據(jù)。

數(shù)值優(yōu)化方法參數(shù)設(shè)置：目標(biāo)質(zhì)量、幾何和氣動(dòng)參數(shù)參考文獻(xiàn)[14]，最大動(dòng)壓約束為1.8×105Pa，最大法向過(guò)載約束為6，最大駐點(diǎn)熱流密度約束為6×105W/m2。攻角最大值為30°，傾側(cè)角最大值為85°，分別計(jì)算最大縱程、縱程為最大縱程的0.9、0.8、0.7、0.6、0.5和0.4倍條件下的最大橫程共13條優(yōu)化彈道，對(duì)應(yīng)終點(diǎn)連接構(gòu)成數(shù)值優(yōu)化方法計(jì)算的可達(dá)區(qū)邊界。

常傾側(cè)角方法參數(shù)設(shè)置：攻角采用給定的標(biāo)稱攻角曲線，分別計(jì)算傾側(cè)角為0°、15°、25°、35°、45°、55°和65°的13條彈道，對(duì)應(yīng)終點(diǎn)連接構(gòu)成常傾側(cè)角方法計(jì)算的可達(dá)區(qū)邊界。

將三種方法計(jì)算的可達(dá)區(qū)邊界進(jìn)行比較，如圖2所示。

圖2 三種方法計(jì)算的可達(dá)區(qū)邊界Fig.2 Footprint boundaries by three algorithms

通過(guò)比較可以看出，本文算法計(jì)算的最大縱程和橫程與優(yōu)化結(jié)果相近(最大縱程誤差約1.9%，最大橫程誤差約3%)，可達(dá)區(qū)邊界的橢圓形近似較為合理，可達(dá)區(qū)橢圓邊界幾乎與優(yōu)化的可達(dá)區(qū)邊界重合，由此本文算法的準(zhǔn)確性和有效性得到驗(yàn)證。本文算法結(jié)果接近了滿足過(guò)程約束的數(shù)值優(yōu)化結(jié)果，這是因?yàn)椋孩僮畲笊璞然璧募僭O(shè)保證其最大射程的近似最優(yōu)性；②平衡滑翔的假設(shè)則保證了所算彈道通常滿足過(guò)程約束。通過(guò)比較本文和數(shù)值優(yōu)化方法計(jì)算的最大橫程對(duì)應(yīng)的傾側(cè)角變化規(guī)律(如圖3所示)可以看出，本文構(gòu)造的傾側(cè)角變化曲線在動(dòng)壓較大的滑翔段基本逼近了數(shù)值優(yōu)化結(jié)果，也側(cè)面驗(yàn)證了本文計(jì)算結(jié)果的近似最優(yōu)性。

圖3 最大橫程對(duì)應(yīng)的傾側(cè)角變化曲線Fig.3 Bank angle curves corresponding maximum crossrange

相比之下，常傾側(cè)角方法計(jì)算的可達(dá)區(qū)邊界雖在面積和形狀上與其他兩種結(jié)果相近，但存在向初始點(diǎn)方向的整體性偏移，這說(shuō)明：①常傾側(cè)角45°的控制律確實(shí)可以近似獲得最大橫程，具有較高的工程應(yīng)用價(jià)值；②攻角采用標(biāo)稱曲線，縱程沒(méi)有優(yōu)化，導(dǎo)致其射程偏小，造成了可達(dá)區(qū)邊界的整體性前移。

在相同軟硬件環(huán)境(Core i7 雙核處理器，MATLAB2016軟件)下，對(duì)三種算法耗時(shí)進(jìn)行測(cè)量，結(jié)果為：數(shù)值優(yōu)化方法約231 s，常傾側(cè)角方法約1.04 s，本文算法0.24 s。可以看出，數(shù)值優(yōu)化方法耗時(shí)最長(zhǎng)，且依賴目標(biāo)的先驗(yàn)信息量較大，不太適用于防御方的目標(biāo)實(shí)時(shí)預(yù)警；常傾側(cè)角方法和本文方法由于不需要優(yōu)化迭代，耗時(shí)較短，但由于常傾側(cè)角方法計(jì)算彈道條數(shù)較多，因此其耗時(shí)略大于本文方法。

本文算法依賴較少的目標(biāo)先驗(yàn)信息，且具有較高的運(yùn)算精度和速度，可對(duì)HGRV類目標(biāo)進(jìn)行準(zhǔn)實(shí)時(shí)連續(xù)的可達(dá)區(qū)預(yù)測(cè)，適用于防御方針對(duì)HGRV類目標(biāo)的實(shí)時(shí)預(yù)警。隨著HGRV飛行速度的逐漸減小，所計(jì)算的可達(dá)區(qū)橢圓邊界精度會(huì)逐漸提高，可達(dá)區(qū)面積逐漸收縮，為目標(biāo)威脅評(píng)估和防御決策提供的信息也更為準(zhǔn)確。另外，可達(dá)區(qū)橢圓邊界計(jì)算的準(zhǔn)確程度也強(qiáng)依賴于前期對(duì)HGRV狀態(tài)參數(shù)和氣動(dòng)參數(shù)估計(jì)的準(zhǔn)確程度。

6 結(jié)論

本文基于平衡滑翔和最優(yōu)化飛行準(zhǔn)則，在獲得有限的HGRV彈道估計(jì)參數(shù)(基于雷達(dá)探測(cè)數(shù)據(jù)通過(guò)實(shí)時(shí)彈道估計(jì)獲得目標(biāo)當(dāng)前時(shí)刻的位置、速度和最大升阻比參數(shù))條件下，提出了一種適用于防御方進(jìn)行HGRV類目標(biāo)預(yù)警的可達(dá)區(qū)邊界快速計(jì)算方法。仿真結(jié)果表明，與傳統(tǒng)方法相比，本文方法具有利用先驗(yàn)信息少、精度較高和運(yùn)算量小的優(yōu)點(diǎn)，可以滿足實(shí)時(shí)性需求，可為針對(duì)HGRV類目標(biāo)的威脅判斷和防御決策提供有效的數(shù)據(jù)支持。