許強強, 葛健全, 楊 濤, 陶 燁, 黃 浩

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

0 引 言

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

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

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

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

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

1 數(shù)學模型

1.1 再入動力學模型

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

(1)

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

(2)

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

1.2 端點約束

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

(3)

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

(4)

1.3 控制變量約束

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

(5)

1.4 過程約束模型

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

(6)

(7)

(8)

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

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

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

Rk∞S/N

(9)

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

理想情況下雷達的信噪比為

(10)

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

由雷達方程可以求出:

(11)

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

令

則有

進一步有

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

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

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

(12)

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

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

(13)

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

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

(14)

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

(15)

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

2 求解策略

2.1 hp-自適應Radau偽譜法

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

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

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

(16)

式中,Φ為非積分項指標;g為積分項指標的被積分項。

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

2.2 逐步增加約束方法

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

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

其具體步驟為:

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

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

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

3 仿真結果及分析

3.1 仿真條件

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

(17)

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

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

表1 多約束條件設置

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

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

Jtime=tf

(18)

(19)

(1) 繞飛情況對比

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

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

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

Jmin time=3.773 1e-18

(20)

Jmin risk=8.646 1e-22

(21)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結 論

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

參考文獻：

[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] 陳小慶,侯中喜,劉建霞. 高超聲速滑翔飛行器再入軌跡多目標多約束優(yōu)化[J]. 國防科技大學學報,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]. 宇航學報,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] 張夢櫻, 唐乾剛, 韓小軍, 等. 復雜約束條件下的再入軌跡迭代求解方法[J]. 兵工學報,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] 丁鷺飛,耿富錄.雷達原理[M].西安:西安電子科技大學出版社, 2003: 128-134.

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

[16] 王海濤,李軍營,梁立威,等.基于hp自適應Radau偽譜法的再入飛行器軌跡優(yōu)化[J].科學技術與工程,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] 劉鶴鳴,丁達理,黃長強,等.基于自適應偽譜法的UCAV低可探測攻擊軌跡規(guī)劃研究[J].系統(tǒng)工程與電子技術,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].宇航學報, 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.