面向突防的多約束滑翔彈道優(yōu)化設(shè)計(jì)

許強(qiáng)強(qiáng), 葛健全, 楊 濤, 陶 燁, 黃 浩

(國防科技大學(xué)航天科學(xué)與工程學(xué)院, 湖南 長沙 410073)

0 引 言

在導(dǎo)彈突防過程中,隱身性能發(fā)揮著重要的作用,提高隱身技術(shù)一方面可從降低自身特征信號入手,另一方面則是通過彈道規(guī)劃,降低雷達(dá)探測概率。因此,采取快速通過雷達(dá)探測區(qū)、反導(dǎo)防御區(qū)的戰(zhàn)術(shù)規(guī)避技術(shù),能夠有效增加高超聲速助推滑翔飛行器打擊目標(biāo)的突然性,提高突防性能。如何在彈道設(shè)計(jì)過程中描述此類路徑(本文定義為輻射型禁飛區(qū))約束則是首要解決問題。

目前對于路徑約束的研究主要側(cè)重于禁飛區(qū)建模方面,Timothy等[1]在平面地球的模型下,建立了無限長圓柱的禁飛區(qū)模型并設(shè)計(jì)了最優(yōu)軌跡;在此基礎(chǔ)上,相關(guān)學(xué)者對該禁飛區(qū)模型下的軌跡優(yōu)化[2-5]、制導(dǎo)[6-9]等方面內(nèi)容進(jìn)行了研究。陳小慶等[10]則是基于球面模型的基礎(chǔ)上利用球面三角形原理建立了無限長圓柱的禁飛區(qū)模型;謝愈等[11]建立了半橢球的禁飛區(qū)模型,張夢櫻[12]基于立體幾何原理在地心系下建立了禁飛區(qū)的半橢球空間幾何數(shù)學(xué)模型。

這樣的模型一方面難以反映飛行器可以進(jìn)入但存在一定的暴露危險(xiǎn)的特征;另一方面絕對的邊界鎖定使得求解更加困難,甚至無法得到最優(yōu)解。同時(shí),在實(shí)際作戰(zhàn)過程中,雷達(dá)探測區(qū)域半徑較大,采用上述模型將對飛行器橫側(cè)向機(jī)動能力提出嚴(yán)峻挑戰(zhàn)。

為更準(zhǔn)確方便地處理輻射型禁飛區(qū)約束,本文根據(jù)將禁飛區(qū)內(nèi)威脅量化的思想[13],將禁飛區(qū)約束轉(zhuǎn)化為與威脅系數(shù)積分相關(guān)的目標(biāo)函數(shù),并通過分析雷達(dá)探測模型,建立輻射型禁飛區(qū)內(nèi)的威脅模型,然后得到輻射型禁飛區(qū)約束對應(yīng)的最優(yōu)控制問題模型。

最后考慮到實(shí)際飛行距離長、空域廣、速域?qū)挼忍攸c(diǎn),飛行器受到嚴(yán)峻的動壓、過載、熱流等過程約束的影響,飛行軌跡將被限制在較為狹窄的范圍。對于此種同時(shí)含有端點(diǎn)約束、過程約束、路徑約束等多約束下的軌跡優(yōu)化問題可以歸結(jié)為強(qiáng)非線性、多階段、多約束的最優(yōu)控制問題,其求解十分困難。偽譜法作為求解復(fù)雜最優(yōu)控制問題的工具,近年來得到迅速發(fā)展。本文基于hp-自適應(yīng)Radau偽譜法、借鑒逐步增加約束方法完成對此軌跡優(yōu)化問題求解,得到優(yōu)化的滑翔段彈道。

1 數(shù)學(xué)模型

1.1 再入動力學(xué)模型

忽略地球的自轉(zhuǎn)和非球形影響,在極坐標(biāo)系下建立高超聲速滑翔飛行器的再入運(yùn)動方程[14]:

(1)

式中,m為飛行器質(zhì)量;g為重力加速度;r=R+h為地心距,r為地球半徑,h為飛行器高度;λ為地理經(jīng)度;φ為地理緯度;V為飛行器速度;θ為當(dāng)?shù)厮俣葍A角;σ為速度偏角(參考方向?yàn)楫?dāng)?shù)卣狈较?;ν為飛行器的傾側(cè)角;L為升力;D為阻力。升力L和阻力D的表達(dá)式分別為

(2)

式中,ρ為空氣密度;ρV2/2為動壓頭;ρ0e-β(r-r0)為指數(shù)大氣模型;Cl為升力系數(shù);Cd為阻力系數(shù);S為飛行器特征面積。

1.2 端點(diǎn)約束

對于初始狀態(tài),通常為固定值;而終端約束通常是由飛行任務(wù)和末制導(dǎo)交接班要求所決定,若要使飛行器到達(dá)指定點(diǎn),則需對位置參數(shù):地心距、地理經(jīng)度、地理緯度進(jìn)行約束,即

(3)

根據(jù)作戰(zhàn)任務(wù)要求,還需要對末端攻擊速度和攻擊角度進(jìn)行限制,即

(4)

1.3 控制變量約束

控制變量取為飛行器的攻角α和傾側(cè)角ν,其約束為

(5)

1.4 過程約束模型

飛行器再入過程約束通常包括:動壓約束、法向過載約束、駐點(diǎn)熱流約束。根據(jù)基本的氣動力熱的理論,結(jié)合相關(guān)的工程估算技術(shù),得到過程約束模型如下:

(6)

(7)

(8)

1.5 威脅模型與性能指標(biāo)函數(shù)

針對輻射型禁飛區(qū)內(nèi)不同的空間位置所對應(yīng)的威脅不同,提出如下假設(shè):任意位置上的威脅可以通過單位時(shí)間內(nèi)該位置上的威脅系數(shù)來量化。相應(yīng)的問題就轉(zhuǎn)化為如何確定禁飛區(qū)內(nèi)任意位置的威脅系數(shù),即禁飛區(qū)的威脅模型。

相關(guān)研究[15]表明,雷達(dá)發(fā)現(xiàn)目標(biāo)的概率與目標(biāo)位置對應(yīng)的信噪相關(guān),假設(shè)在雷達(dá)覆蓋區(qū)內(nèi),任意位置的威脅系數(shù)與該位置的信噪比成正比,即:

Rk∞S/N

(9)

式中,Rk為威脅系數(shù);S/N為雷達(dá)的信噪比。

理想情況下雷達(dá)的信噪比為

(10)

式中,Ps是無氣象干擾時(shí)的接收信號功率;N是雷達(dá)噪聲功率,視為常數(shù)。

由雷達(dá)方程可以求出:

(11)

式中,Pt,G,δ,λ分別為雷達(dá)的發(fā)射機(jī)功率、天線的增益、目標(biāo)的雷達(dá)截面積以及工作波長,在目標(biāo)飛行器和雷達(dá)都確定的情況下都可以視為常數(shù);Rd為理想情況下的雷達(dá)作用距離,即雷達(dá)與目標(biāo)飛行器的距離。

令

則有

進(jìn)一步有

其中當(dāng)雷達(dá)和目標(biāo)飛行器都確定的情況下,K和N均可視為常數(shù),理想情況下信噪比與距離的關(guān)系如圖1所示。

圖1 理想情況下信噪比與作用距離的關(guān)系圖Fig.1 Signal to noise ratio vs. distance in ideal conditions

又威脅系數(shù)Rk∞S/N,所以

(12)

式中,C為比例系數(shù),可根據(jù)實(shí)際情況由用戶確定。

為了便于計(jì)算,定義威脅系數(shù)因子:

(13)

理想情況下,在雷達(dá)和目標(biāo)飛行器都確定的情況下,Cr為常數(shù),由雷達(dá)系統(tǒng)確定。

得到威脅系數(shù)Rk表達(dá)式如下:

(14)

這里,為了計(jì)算的方便,假設(shè)雷達(dá)的探測距離為無窮遠(yuǎn),當(dāng)Rd很大,則Rk很小,其對積分結(jié)果的影響很小,故可以對Rk進(jìn)行全彈道積分,由此得到基于威脅系數(shù)積分的彈道優(yōu)化性能指標(biāo)函數(shù)如下:

(15)

從該函數(shù)可以看出,全程威脅積分越小,則需飛行器與雷達(dá)的距離越遠(yuǎn),相應(yīng)地會增加飛行距離;而飛行距離的增加,使得飛行時(shí)間延長,全程威脅積分增加。因此,存在一個(gè)飛行時(shí)間和作用距離的組合使得全程威脅積分最小。顯然,該性能指標(biāo)函數(shù)成功地將飛行時(shí)間和繞飛距離兩個(gè)指標(biāo)綜合到了一起,為輻射型禁飛區(qū)的處理提供了便利。

2 求解策略

2.1 hp-自適應(yīng)Radau偽譜法

近年來,由于偽譜法具有收斂半徑大、對初值不敏感、收斂速度快等特點(diǎn),被廣泛應(yīng)用于飛行器軌跡優(yōu)領(lǐng)域中,并取得了一定的研究成果[16]。hp-自適應(yīng)Radau偽譜法結(jié)合了優(yōu)化網(wǎng)格數(shù)目(h方法的計(jì)算稀疏性)和增階(p方法的快速收斂性)兩方面的優(yōu)勢,能夠較好地處理具有多約束條件下的最優(yōu)控制問題[17]。Radau偽譜法的基本原理是:將每一個(gè)網(wǎng)格區(qū)間上的狀態(tài)變量和控制變量在眾多LGR(Legendre-Gauss-Radau)點(diǎn)上進(jìn)行離散,并將離散點(diǎn)作為節(jié)點(diǎn),進(jìn)而創(chuàng)建Lagrange插值多項(xiàng)式,以擬合近似狀態(tài)變量和控制變量。通過對全局插值多項(xiàng)式得到的狀態(tài)變量進(jìn)行求導(dǎo),來近似狀態(tài)變量對時(shí)間的導(dǎo)數(shù),從而可以將系統(tǒng)的動力學(xué)微分方程約束轉(zhuǎn)化為代數(shù)約束,性能指標(biāo)中的積分則由Radau積分進(jìn)行計(jì)算。由于全局當(dāng)前網(wǎng)格區(qū)間的終端狀態(tài)為下一網(wǎng)格區(qū)間的初始狀態(tài),由此避免了全局Radau偽譜法中終端狀態(tài)的積分過程。通過上述一系列步驟,可以將原有的最優(yōu)控制問題變?yōu)榉蔷€性規(guī)劃問題(nonlinear programming,NLP)。

Radau偽譜法的求解具體步驟可參考相關(guān)文獻(xiàn)[18-19],本文不再詳述。其中將Bolza型性能指標(biāo)函數(shù)中的積分項(xiàng)用LGR積分近似,得到近似性能指標(biāo)函數(shù)如下:

J=Φ(X0,t0,Xf,tf)+

(16)

式中,Φ為非積分項(xiàng)指標(biāo);g為積分項(xiàng)指標(biāo)的被積分項(xiàng)。

對于本文,性能指標(biāo)函數(shù)如式(15)所示。

2.2 逐步增加約束方法

考慮多約束條件下飛行器軌跡優(yōu)化問題比較復(fù)雜,直接采用偽譜法求解難以得到優(yōu)化結(jié)果。故借鑒于文獻(xiàn)[12],采用逐步增加約束的方法,將具有較少約束問題的優(yōu)化結(jié)果作為下一步的優(yōu)化初值進(jìn)行迭代,流程如圖2所示。

圖2 逐步增加約束方法流程圖Fig.2 Flow chart of adding constraints

其具體步驟為:

步驟1基于前面建立的最優(yōu)控制問題的模型,添加必要的端點(diǎn)約束并配置好數(shù)值算法,以最短時(shí)間為優(yōu)化目標(biāo),得到滿足端點(diǎn)約束的最優(yōu)彈道,驗(yàn)證目標(biāo)可達(dá)性;

步驟2在步驟1的模型基礎(chǔ)上添加飛行過程約束,以步驟1的優(yōu)化結(jié)果作為求解初值,以最短時(shí)間為優(yōu)化目標(biāo),保證在過程約束下目標(biāo)可達(dá)性;

步驟3在步驟2的基礎(chǔ)上,添加輻射型禁飛區(qū)約束,將性能指標(biāo)函數(shù)改成威脅最小,以步驟2中最短時(shí)間彈道為初值,求解繞過禁飛區(qū)時(shí)的威脅最小彈道。

3 仿真結(jié)果及分析

3.1 仿真條件

取飛行器質(zhì)量為1 500 kg,參考面積1 m2。嚴(yán)格來講,氣動系數(shù)通常與飛行高度、馬赫數(shù)和攻角三者相關(guān),可表示為三者的函數(shù)。但在高超聲速條件下,飛行器的氣動系數(shù)遵循馬赫數(shù)無關(guān)原理[20],可以認(rèn)為氣動系數(shù)僅隨攻角變化,本文中取飛行器的氣動系數(shù)與攻角的關(guān)系為

(17)

式中,Cl,Cd分別為飛行器的升力系數(shù)和阻力系數(shù);α為飛行器攻角。

具體仿真參數(shù)如表1所示。同時(shí),由于Cr=1,計(jì)算過程中,為了降低計(jì)算難度、提高計(jì)算效率,本文采用無量綱化的方法,取長度的單位長度為地球半徑,取質(zhì)量的單位質(zhì)量為飛行器質(zhì)量。

表1 多約束條件設(shè)置

3.2 優(yōu)化結(jié)果及分析

本文仿真選擇最短時(shí)間彈道和最小威脅彈道作對比。優(yōu)化目標(biāo)分別為

Jtime=tf

(18)

(19)

(1) 繞飛情況對比

最短時(shí)間彈道和最小威脅彈道的結(jié)果如圖3所示,實(shí)線表示的是最短時(shí)間彈道,虛線表示的是最小威脅彈道,可以很明顯的看出飛行器成功繞過禁飛區(qū)中心。

圖3 最短時(shí)間彈道和威脅最小彈道對比圖Fig.3 Comparison diagram of minimum time trajectory and minimum threat trajectory

對時(shí)間最短彈道和威脅最小彈道的威脅系數(shù)積分,結(jié)果分別為

Jmin time=3.773 1e-18

(20)

Jmin risk=8.646 1e-22

(21)

從積分結(jié)果對比可以看到威脅最小彈道的威脅系數(shù)積分與最短時(shí)間彈道的相差10-4量級,從這一方面也說明了這種輻射型禁飛區(qū)的處理方式是正確有效的。

通過簡單的理論分析可以得到,對于固定中心的輻射型禁飛區(qū),由于威脅系數(shù)與距離的關(guān)系,若使全程的威脅積分最小,可在兩個(gè)方面做出變化:一是增大飛行器與禁飛區(qū)中心的徑向距離,二是考慮到時(shí)間對積分的影響,當(dāng)威脅系數(shù)相對較大時(shí),提高飛行器速度,縮短通過高威脅區(qū)域時(shí)間。

從圖3(b)可以看出,飛行器比較靠近禁飛區(qū)時(shí)的經(jīng)度范圍為35°～45°,兩條彈道對應(yīng)的飛行時(shí)間范圍為500～900 s。

下面將結(jié)合圖中兩彈道信息的對比分析上述縮小全程威脅積分的2個(gè)方面分別的體現(xiàn):

①從圖3(a)的平面軌跡對比中可以很明顯的看出,在接近禁飛區(qū)中心時(shí),最小威脅彈道與中心的水平距離明顯增大,相對減小了威脅系數(shù);

②從圖3(c)的高度時(shí)間變化關(guān)系對比中可以看出,在500～900 s的范圍內(nèi),最小威脅彈道的高度總體要比最短時(shí)間彈道的高,但是差別無水平距離明顯;

③從圖3(d)的速度時(shí)間變化關(guān)系對比,顯然地最小威脅彈道的速度要高于最短時(shí)間彈道,對應(yīng)的就是如果在覆蓋區(qū)內(nèi)飛過同樣的距離,最小威脅彈道用時(shí)更短,威脅積分也更小。

通過上述分析,最小威脅彈道的仿真結(jié)果所體現(xiàn)的繞飛途徑與理論分析相吻合,這也進(jìn)一步驗(yàn)證了本文中建立的輻射型禁飛區(qū)的威脅模型和將輻射型禁飛區(qū)約束轉(zhuǎn)換為與威脅積分相關(guān)的性能指標(biāo)函數(shù)的處理方法的正確性。

(2) 約束滿足度和可行性分析

最小威脅彈道的過程約束和控制變量隨時(shí)間的變化關(guān)系如圖4所示。

圖4 最小威脅彈道的過程約束和控制變量變化曲線Fig.4 Curves of process constraints and control variables in minimum threat trajectory

關(guān)于過程約束,由圖4(a)～圖4(c)可知,法向過載、駐點(diǎn)熱流和動壓均滿足給定的約束范圍。

關(guān)于控制變量,攻角大部分時(shí)間穩(wěn)定在4°附近,僅在結(jié)束時(shí)變化較大,這也是符合飛行器滑翔段保持小攻角飛行的特征;傾側(cè)角在全程為負(fù)值,在繞飛階段達(dá)到極值,這也符合飛行軌跡左傾且在繞飛階段曲率較大的特征。

關(guān)于終端約束,最小威脅彈道能準(zhǔn)確到達(dá)終點(diǎn)(60°E,35°N),且Hf=20 km,同時(shí)速度也滿足給定的終端速度范圍。

綜上,最小威脅彈道能夠滿足所給定的過程約束,同時(shí)其優(yōu)化得到的控制變量的變化在給定約束范圍之間,故最小威脅彈道是可行的。

4 結(jié) 論

本文面向助推滑翔飛行器突防,對多約束條件下的滑翔段彈道進(jìn)行了優(yōu)化設(shè)計(jì)。針對傳統(tǒng)路徑約束建模的不足,建立了輻射型禁飛區(qū)模型。采用了逐步增加約束的方法,利用hp-自適應(yīng)Radau偽譜法對滑翔段彈道進(jìn)行了優(yōu)化設(shè)計(jì),成功得到了滿足相關(guān)約束且使全程威脅積分最小的優(yōu)化彈道,達(dá)到了快速突防的目的。并且結(jié)果彈道反映的信息與理論分析一致,也證明了本文中建立的威脅模型和針對輻射型禁飛區(qū)時(shí)處理方法的正確性,為面向突防的滑翔段彈道優(yōu)化設(shè)計(jì)提供了有效的方法。

參考文獻(xiàn)：

[1] JORRIS T R, COBB R G. Three-dimensional trajectory optimization satisfying waypoint and no-fly zone constraints[J].Journal of Guidance, Control, and Dynamics, 2009, 33(2): 551-572.

[2] WANG L, XING Q H, MAO Y F. Reentry trajectory rapid optimization for hypersonic vehicle satisfying waypoint and no-fly zone constraints[J]. Journal of Systems Engineering and Electronics, 2015, 26(6): 1277-1290.

[3] ZHAO D J, SONG Z Y. Reentry trajectory optimization with waypoint and no-fly zone constraints using multiphase convex programming[J].Acta Astronautica, 2017, 137: 60-69.

[4] MAO Y F, ZHANG D L, WANG L. Reentry trajectory optimization for hypersonic vehicle based on improved Gauss pseudospectral method[J]. Soft Computing, 2016, 21(16): 1-10.

[5] LIU X F, SHEN Z J, LU P. Entry trajectory optimization by second-order cone programming[J]. Journal of Guidance Control & Dynamics, 2015, 39(2): 1-15.

[6] LIANG Z X, LIU S Y, LI Q D, et al. Lateral entry guidance with no-fly zone constraint[J].Aerospace Science & Technology,2017,60: 39-47.

[7] ZHANG D, LIU L, WANG Y J. On-line reentry guidance algorithm with both path and no-fly zone constraints[J]. Acta Astronautica, 2015,117: 243-253.

[8] LIANG Z X, REN Z, BAI C. Lateral reentry guidance for maneuver glide vehicles with geographic constraints[C]∥Proc.of the Control Conference, 2013: 5187-5192.

[9] GUO J, WU X Z, TANG S J. Autonomous gliding entry guidance with geographic constraints[J]. Chinese Journal of Aeronautics, 2015, 28(5): 1343-1354.

[10] 陳小慶,侯中喜,劉建霞. 高超聲速滑翔飛行器再入軌跡多目標(biāo)多約束優(yōu)化[J]. 國防科技大學(xué)學(xué)報(bào),2009,31(6):77-83.

CHEN X Q, HOU Z X, LIU J X. Multi-objective optimization of reentry trajectory for hypersonic glide vehicle with multi-constraints[J].Journal of National University of Defense Technology, 2009, 31(6): 77-88.

[11] 謝愈, 劉魯華, 湯國建. 多約束條件下高超聲速度滑翔飛行器軌跡優(yōu)化[J]. 宇航學(xué)報(bào),2011,32(12):2499-2504.

XIE Y, LIU L H, TANG G J. Trajectory optimization for hypersonic glide vehicle with multi-constraints[J]. Journal of Astronatics, 2011, 32(12): 2499-2504.

[12] 張夢櫻, 唐乾剛, 韓小軍, 等. 復(fù)雜約束條件下的再入軌跡迭代求解方法[J]. 兵工學(xué)報(bào),2015, 36(6): 1015-1023.

ZHANG M Y, TANG Q G, HAN X J, et al. Iterative method to solving re-entry trajectory optimization with complex constraints[J].Acta armamentarii, 2015, 36(6):1015-1023.

[13] VIAN J L, MOORE J R. Trajectory optimization with risk minimization for military aircraft[J]. Guidance, 1989, 12(3): 311-317.

[14] VINH N X. Optimal trajectories in atmospheric flight[M]. New York: Elsevier Scientific Publishing Company, 1981.

[15] 丁鷺飛,耿富錄.雷達(dá)原理[M].西安:西安電子科技大學(xué)出版社, 2003: 128-134.

DING L F, GENG F L. Principle of radar[M].Xi’an: Xidian University Press, 2003: 128-134.

[16] 王海濤,李軍營,梁立威,等.基于hp自適應(yīng)Radau偽譜法的再入飛行器軌跡優(yōu)化[J].科學(xué)技術(shù)與工程,2015,15(2):165-171.

WANG H T, LI J Y, LIANG L W, et al. Track optimizing for reentry vehicle based on hp- adaptive radau pseudospectural method[J]. Science Technology and Engineering,2015,15(2): 165-171.

[17] 劉鶴鳴,丁達(dá)理,黃長強(qiáng),等.基于自適應(yīng)偽譜法的UCAV低可探測攻擊軌跡規(guī)劃研究[J].系統(tǒng)工程與電子技術(shù),2013,35(1): 78-84.

LIU H M, DING D L, HUANG C Q, et al. UCAV low observable attacking trajectory planning based on adaptive pseudospectral method[J].Systems Engineering and Electronics,2013,35(1):78-84.

[18] CHRISTOPHER L D, WILLIAM W H, ANIL V R. An hp-adaptive pseudospectral method for solving optimal problems[J].Optimal Control Applications and Methods,2011,32:476-502.

[19] CAMILA C F, ANIL V R. Direct trajectory optimization and costate estimation of state inequality path-constrained optimal control problems using a radau pseudospectral method[C]∥Proc.of the AIAA Guidance,Navigation,and Control Confe-rence, 2012.

[20] 雍恩米,唐國金,陳磊.基于Gauss偽譜方法的高超聲速飛行器再入軌跡快速優(yōu)化[J].宇航學(xué)報(bào), 2008, 29(6): 1766-1772.

YONG E M, TANG G J, CHEN L. Rapid trajectory optimization for hypersonic reentry vehicle via gauss psedospectral method[J].Journal of Astronautics,2008,29(6):1766-1772.