單脈沖最省能量攔截軌道迭代算法

孟少飛, 劉新學(xué), 傅 丹, 楊 其

(1.火箭軍工程大學(xué)初級(jí)指揮學(xué)院, 陜西 西安 710025; 2. 中國(guó)航天科技集團(tuán)公司第一研究院, 北京 100071)

?

單脈沖最省能量攔截軌道迭代算法

孟少飛1, 劉新學(xué)1, 傅 丹2, 楊 其1

(1.火箭軍工程大學(xué)初級(jí)指揮學(xué)院, 陜西 西安 710025; 2. 中國(guó)航天科技集團(tuán)公司第一研究院, 北京 100071)

最優(yōu)攔截軌道解算是航天器在軌攔截問(wèn)題中的重要內(nèi)容,是天基武器研發(fā)中需要解決的核心技術(shù)難題之一。由于軌道解算方程和攔截方程是隱式非線性方程,因此很難用解析的方法進(jìn)行分析和計(jì)算。針對(duì)非固定攔截點(diǎn)條件下的攔截軌道優(yōu)化問(wèn)題,分析討論了攔截軌道參數(shù)和攔截器變軌點(diǎn)特征速度隨攔截時(shí)間和攔截真近點(diǎn)角的變化規(guī)律。在此基礎(chǔ)上,研究了由攔截軌道的真近點(diǎn)角反向求解攔截初始參數(shù)問(wèn)題,并設(shè)計(jì)了能量最省攔截軌道的快速迭代算法。最后通過(guò)仿真,驗(yàn)證了算法的可行性和快速性。

空間攔截; 特征速度; 快速迭代

0 引 言

隨著現(xiàn)代空間技術(shù)的發(fā)展,奪取制天權(quán)并保持空間優(yōu)勢(shì)已經(jīng)成為世界各國(guó)的重要軍事發(fā)展戰(zhàn)略。同時(shí),空間對(duì)抗技術(shù)與手段隨著空間技術(shù)的發(fā)展進(jìn)步而變得越來(lái)越復(fù)雜與多樣,其約束性越來(lái)越多,偶然性也越來(lái)越大?？臻g攔截作為空間對(duì)抗的重要內(nèi)容,其首要的約束條件即是攔截器所攜帶的能量的約束[1-2],除此之外還有其他戰(zhàn)術(shù)條件上的約束,包括攔截時(shí)間[3-4]、攔截窗口[5-6]、攔截位置[7]等約束條件。隨著多脈沖變軌技術(shù)和多圈攔截交會(huì)技術(shù)的發(fā)展,雖然多脈沖攔截和多圈攔截所需的能量往往小于單脈沖單圈攔截所需的能量,但是考慮到天基武器的特殊作戰(zhàn)需求以及“發(fā)出后不管”的高效作戰(zhàn)模式,單脈沖攔截相比多脈沖攔截能減小作戰(zhàn)所需的后續(xù)軌道保障成本和可靠性;而相比于多圈攔截模式,單圈攔截能較少攔截軌道和目標(biāo)軌道的交叉,減小暴露作戰(zhàn)意圖和軌道被跟蹤探測(cè)的風(fēng)險(xiǎn)。

由于最優(yōu)攔截軌道的解算問(wèn)題是一個(gè)具有多個(gè)局部最優(yōu)解的非線性優(yōu)化問(wèn)題,如何避免算法陷入局部最優(yōu)成了首先需要考慮的問(wèn)題。文獻(xiàn)[8-9]研究了起點(diǎn)和終點(diǎn)固定而飛行時(shí)間不定情況下所需能量隨飛行時(shí)間的變化規(guī)律和極小能量軌道;文獻(xiàn)[10-11]采用了全局搜索的方法實(shí)現(xiàn)了對(duì)單脈沖最省能量攔截軌道的計(jì)算;文獻(xiàn)[12-17]以能量最省、攔截時(shí)間最短、脫靶量最小為綜合指標(biāo),采用遺傳算法等組合優(yōu)化算法對(duì)單脈沖攔截軌道的優(yōu)化進(jìn)行了研究。而本文考慮在軌攔截中攔截起點(diǎn)固定而攔截終點(diǎn)不固定條件下的單脈沖單圈最省能量攔截問(wèn)題,從轉(zhuǎn)移速度增量的變化規(guī)律角度入手設(shè)計(jì)了能量最省攔截軌道的數(shù)值迭代算法,既提高了算法的快速性與可靠性,又避免了陷入局部最優(yōu)的問(wèn)題。

1 非固定終端攔截特征速度變化規(guī)律

轉(zhuǎn)移軌道示意圖如圖1所示,航天器先后通過(guò)固定端點(diǎn)P1、P2時(shí)的位置矢量為r1、r2,速度矢量為v1、v2,c為連接P1、P2的弦長(zhǎng),當(dāng)其軌道為橢圓時(shí),所需要的時(shí)間Δt和轉(zhuǎn)移軌道半長(zhǎng)軸a之間的關(guān)系[18-19]為

(1)

這就是Lagrange轉(zhuǎn)移時(shí)間方程,其中α、β為拉格朗日參數(shù),且有

(2)

(3)

(4)

圖1 轉(zhuǎn)移軌道示意圖Fig.1 Transfer orbit sketch map

1.1 特征速度關(guān)于軌道半長(zhǎng)軸的導(dǎo)數(shù)

設(shè)攔截器變軌點(diǎn)速度增量的大小為Δv,則有

(5)

式中,v+和v-為攔截器變軌前和變軌后的速度；φ為它們的夾角?？梢钥闯?當(dāng)變軌時(shí)刻和攔截時(shí)刻都不確定時(shí),Δv為變量v+、v-和φ的函數(shù);而當(dāng)變軌時(shí)刻確定而攔截時(shí)刻不確定時(shí),Δv僅為變量v+和φ的函數(shù)。

設(shè)a為變軌后攔截軌道的半長(zhǎng)軸,將式(5)兩邊分別對(duì)a求導(dǎo)[20]有

(6)

其中

(7)

(8)

由橢圓軌道的速度公式得

(9)

對(duì)式(1)～式(9)兩邊求導(dǎo)可得

(10)

此時(shí),式(6)右邊4個(gè)因式中,還有dφ/da沒(méi)有解決。為解決這一因式,需要確定φ和變軌后的方向角γ+的關(guān)系,如圖2所示。

圖2 變軌前后速度變化示意圖Fig.2 Change of velocity before and after orbital transfer

圖中,γ-和γ+為變軌前后的飛行方向角,且有γ-<90°,γ+>90°。過(guò)A點(diǎn)作AB垂直于變軌點(diǎn)位置矢量r于B點(diǎn),過(guò)B點(diǎn)在攔截軌道平面內(nèi)作垂直v-的直線交v+于C點(diǎn)。通過(guò)三角函數(shù)關(guān)系可以得到

cosφ=cosγ+cosγ-+sinγ+sinγ-cosθ

(11)

將式(11)兩邊對(duì)a求導(dǎo)可得

(12)

式(11)和式(12)中的γ+、γ-和θ可以由以式(13)求得:

(13)

(14)

(15)

(16)

又由式(13)可得

(17)

(18)

同時(shí),對(duì)式(1)求導(dǎo)可得

(19)

(20)

(21)

至此,可以求出特征速度Δv對(duì)攔截軌道半長(zhǎng)軸a的導(dǎo)數(shù)。

1.2 特征速度的變化規(guī)律

當(dāng)攔截時(shí)間較長(zhǎng)時(shí),可以分為多圈攔截和單圈攔截,多圈攔截是指攔截器變軌后經(jīng)過(guò)多個(gè)轉(zhuǎn)移軌道的軌道周期后攔截指定目標(biāo),單圈攔截則是指攔截器變軌后在一個(gè)軌道周期內(nèi)攔截指定目標(biāo)。顯然,多圈攔截的特征速度往往小于單圈攔截的特征速度。但是多圈攔截的缺點(diǎn)是攔截器多次經(jīng)過(guò)攔截點(diǎn),這樣則容易暴露其攔截意圖,且容易被目標(biāo)航天器跟蹤與規(guī)避。當(dāng)攔截的時(shí)間控制在合適的范圍內(nèi)時(shí),采用單圈攔截需要的最小速度增量和多圈攔截最小速度增量接近,且戰(zhàn)術(shù)效果優(yōu)于多圈攔截。僅考慮單圈攔截的情況下,攔截時(shí)間t和變軌特征速度Δv隨攔截軌道半長(zhǎng)軸a的變化規(guī)律如圖3所示。

圖3 t和Δv隨a的變化規(guī)律Fig.3 Rule of t and Δv with the change of a

由圖3可以看出,當(dāng)攔截時(shí)間較小時(shí),所需要的速度增量較大,如圖中AB段曲線所示。由A點(diǎn)以后的一個(gè)時(shí)間周期內(nèi),攔截特征速度有兩個(gè)極小值,而符合條件的全局極小值只有一個(gè)C點(diǎn),它所對(duì)應(yīng)的攔截時(shí)間為D?？梢钥闯?攔截軌道優(yōu)化問(wèn)題是具有多個(gè)局部最優(yōu)解的非線性優(yōu)化問(wèn)題,在不同的區(qū)間內(nèi)存在不同的局部最優(yōu)解,且全局最優(yōu)解是唯一存在的,而求解最優(yōu)攔截軌道問(wèn)題的本質(zhì)則是求解全局最優(yōu)解。

但在圖3中,攔截軌道特征速度在各區(qū)間的變化是非單調(diào)的,且攔截軌道半長(zhǎng)軸和攔截軌道特征速度之間并非單映射的關(guān)系,無(wú)法用數(shù)值迭代的方法求解局部最優(yōu)攔截軌道。為了方便問(wèn)題的解決,考慮用轉(zhuǎn)移軌道的真近點(diǎn)角Δθ取代轉(zhuǎn)移軌道半長(zhǎng)軸a,討論特征速度Δv的變化規(guī)律。

(22)

(23)

圖4 Δθ的變化范圍Fig.4 Variation range of Δθ

圖5 Δv隨Δθ的變化規(guī)律Fig.5 Rule of Δv with the change of Δθ

可以看出:雖然Δv的變化規(guī)律與變軌時(shí)刻的具體條件有關(guān),但由于目標(biāo)航天器的周期性運(yùn)動(dòng),導(dǎo)致Δv在不同的周期內(nèi),具有相似的變化規(guī)律;在目標(biāo)的一個(gè)運(yùn)動(dòng)周期內(nèi),當(dāng)攔截軌道與原軌道方向一致時(shí),Δv存在唯一最小值Δvmin,且以Δvmin為分界,其兩側(cè)的變化分別都是單調(diào)的。

2 最省能量攔截軌道快速迭代算法

2.1 迭代算法的總體思想和基本模型

(Δθ0,Δv0),…,(Δθn-1,Δvn-1),(Δθn,Δvn)

則有

(24)

設(shè)該區(qū)間內(nèi)的局部最優(yōu)解為F,值為(Δθf(wàn),Δvf),則有

(25)

(26)

當(dāng)計(jì)算精度滿足要求時(shí),迭代停止:

|Δvn-Δvn-1|≤δ

(27)

由式(26)可以看出,迭代條件是至少需要已知前面兩個(gè)點(diǎn),因此迭代初始可以手動(dòng)選擇一個(gè)點(diǎn)進(jìn)行計(jì)算,作為迭代初始條件。當(dāng)獲得多個(gè)時(shí)間周期內(nèi)的多個(gè)局部最優(yōu)解后(當(dāng)N≠0時(shí),取ΔtN=Δt0+TT·N即可),通過(guò)比較每個(gè)局部最優(yōu)解得能量成本和時(shí)間成本,結(jié)合規(guī)定的攔截時(shí)間窗口,從中選擇最合適的解作為全局最優(yōu)解。一般計(jì)算2～3個(gè)周期即可找到全局最優(yōu)解。

迭代算法的流程圖如圖6所示。

圖6 迭代算法流程圖Fig.6 Iteration algorithm flow chart

2.2 最省能量攔截軌道參數(shù)解算

當(dāng)?shù)惴ù_定后,還需要解決如何由攔截軌道的真近點(diǎn)角Δθ計(jì)算攔截器變軌所需要的特征速度Δv。設(shè)m為攔截軌道的方向,1表示順行,-1表示逆行;設(shè)變量Δθ*、m*為攔截軌道所對(duì)應(yīng)的最便捷路徑參數(shù),并與原始參數(shù)Δθ*、m*的關(guān)系為

(28)

(29)

(30)

為了便于求解方程(30),設(shè)

(31)

又設(shè)

(32)

則方程(30)的解可以表示為

(33)

設(shè)方程(30)的兩個(gè)解分別為n1、n2,則需要排除掉一個(gè)不符合實(shí)際情況的解。具體方法是計(jì)算

(34)

(35)

若m*=m′,則對(duì)應(yīng)的ni為所求解。若m*≠m′,則對(duì)應(yīng)的ni需要篩除。當(dāng)n確定后,計(jì)算

(36)

(37)

設(shè)ΔθT為目標(biāo)飛行過(guò)的真近點(diǎn)角度,則有

(38)

3 仿真算例

設(shè)初始時(shí)刻攔截器和目標(biāo)的軌道根數(shù)如表1所示。攔截器在t=0時(shí)刻開(kāi)始機(jī)動(dòng)變軌,攔截時(shí)間約束為Δt≤4 h。通過(guò)仿真,可以得到變軌點(diǎn)速度增量隨攔截器轉(zhuǎn)移真近點(diǎn)角的變化關(guān)系如圖7所示,迭代所得的局部最優(yōu)解如表2中所示。

表1 目標(biāo)和攔截器初始軌道根數(shù)

圖7 仿真算例中Δv隨Δθ的變化規(guī)律Fig.7 The rule of Δv with the change of Δθ in the simulation example

N≤Δθ≤π-Δθ/(°)Δv/(km/s)Δt/hπ+≤Δθ≤2π-Δθ/(°)Δv/(km/s)Δt/h01674．141．773335．520．64511696．174．142911．432．71

通過(guò)對(duì)表2中各個(gè)周期中不同區(qū)間段的局部最優(yōu)解比較可得,當(dāng)攔截軌道真近點(diǎn)角Δθ為291°時(shí),攔截器變軌所需要的特征速度Δv達(dá)到全局最小值1.43 km/s,與之對(duì)應(yīng)的攔截時(shí)間Δt為2.7 h,所需的變軌特征速度矢量為

Δvmin=

對(duì)比圖7中的特征速度變化曲線可以看出,迭代算法求得的最優(yōu)解為全局最優(yōu)解,這驗(yàn)證了迭代算法的正確性。在此條件下的最優(yōu)攔截軌道如圖8所示。

圖8 攔截軌道三維仿真圖Fig.8 The three-dimensional emulational interception orbit

迭代計(jì)算過(guò)程中的迭代次數(shù)和所用時(shí)間如表3所示,可以看出各區(qū)間迭代的平均次數(shù)在30次左右,這對(duì)計(jì)算機(jī)來(lái)說(shuō)屬于較小的計(jì)算量,所用的最大迭代時(shí)間為0.031 s,很好地滿足了天基武器快速性反應(yīng)的要求,同時(shí)也說(shuō)明此迭代方法的快速性。

表3 不同區(qū)間最省能量軌道迭代次數(shù)與時(shí)間

Table 3 The complexity and time of iteration for the optimal interception orbit

N≤Δθ≤π-迭代次數(shù)時(shí)間/sπ+≤Δθ≤2π-迭代次數(shù)時(shí)間/s0300．031310．0311320．016300．025

4 結(jié) 論

文章探討了在一般天基攔截問(wèn)題中,攔截終點(diǎn)不固定條件下的最省能量攔截軌道。首先仿真研究了變軌點(diǎn)速度增量隨攔截器飛過(guò)的真近點(diǎn)角的變化曲線并探討了其共性特點(diǎn);其次在此基礎(chǔ)上研究設(shè)計(jì)了最省能量攔截軌道的分段迭代方法,并進(jìn)行了仿真驗(yàn)證,證明了該方法的有效性。此算法實(shí)現(xiàn)起來(lái)較為簡(jiǎn)單,迭代速度也較快,能滿足星載計(jì)算機(jī)的實(shí)時(shí)計(jì)算能力要求。同時(shí),通過(guò)對(duì)起點(diǎn)固定而終點(diǎn)不固定的最省能量攔截問(wèn)題的研究,能為起點(diǎn)和終點(diǎn)都不固定的無(wú)約束最省能量在軌攔截和交會(huì)問(wèn)題提供一定的借鑒意義和理論支撐。

[1] Henzeh L, Donghoom K, James T. Solution for nonlinear three-dimensional intercept problem with minimum energy[J].MathematicalProblemsinEngineering, 2013: 435725.

[2] Chen X G, Wu Z L, Gao Y. Design and optimization of space-based interception orbit based on energy[J].ComputerSimulation,2012,29(5):59-63.(陳小剛,吳志林,高宇.基于能量的天基攔截軌道設(shè)計(jì)與優(yōu)化[J].計(jì)算機(jī)仿真,2012,29(5):59-63.)

[3] Bai H B, Feng S X, Zhu L P. Simulation of fixed-time orbit interception in space warfare[J].AerospaceControl, 2006, 24(4): 62-66.

[4] Gao H, Zhu Z X, Liu J. An effective minimum-time orbital transfer design under continuous thrust[J].JournalofNorthwesternPolytechnicalUniversity,2012, 30(2): 187-191.(高懷,朱戰(zhàn)霞,劉劍. 基于遺傳算法的連續(xù)推力最短時(shí)間轉(zhuǎn)移軌道設(shè)計(jì)[J]. 西北工業(yè)大學(xué)學(xué)報(bào), 2012, 30(2): 187-191.)

[5] Leeghim H. Spacecraft intercept using minimum control energy and wait time[J].CelestialMechanicsandDynamicalAstronomy, 2013, 115(1): 1-19.

[6] Liu H, Lei T, Li S L, et al. Optimal design for impulsive orbit maneuver based on parameters optimization thinking[J].AerospaceShanghai, 2010(1)：20-25.(劉恒,雷濤,李順利,等.基于參數(shù)優(yōu)化的脈沖軌道設(shè)計(jì)[J].上海航天,2010(1)：20-25.)

[7] Wei P T, Lei G, Guo H N, et al. Orbit interception at virtual crossing point based on Lambert transfer[J].FlightDynamics, 2012, 30(6):556-559.(魏鵬濤,雷剛,郭洪娜,等.基于Lambert轉(zhuǎn)移的虛交點(diǎn)軌道攔截優(yōu)化[J].飛行力學(xué),2012,30(6):556-559.)

[8] Lu T, Zhu J. A genetic algorithm for finding a path subject to two constraints[J].AppliedSoftComputing, 2013, 13(2):891-898.

[9] Ulybyshev Y. Continuous thrust orbit transfer optimization using large-scale linear programming[J].JournalofGuidanceControlandDynamics, 2007, 30(2): 427-436.

[10] Chao T, Wang S Y, Wang Z C, et al. Near space vehicle trajectory optimization approach based on hybrid SVM and GA algorithm[J].JournalofAstronautics, 2012, 33(2): 183-188.(晁濤, 王松艷,王子才,等. 基于組合優(yōu)化算法的臨近空間飛行器軌跡優(yōu)化[J]. 宇航學(xué)報(bào),2012,33(2): 183-188.)

[11] Yan Y H, Guo D C. Indirect optimization for finite-thrust orbital interception problem[J].AppliedMechanicsandMaterials, 2013, 313/314: 1051-1054.

[12] Zheng Q, Sha J X, Shu H, et al. A variant constrained genetic algorithm for solving conditional nonlinear optimal perturbations[J].AdvancesinAtmosphericScience,2014,31(1):219-229.

[13] Leung C S, Lam P M, Tsang P W M. A graphics processing unit accelerated genetic algorithm for affine invariant matching for broken contours[J].JournalofSignalProcessingSystemsforSignalImageandVideoTechndogy,2012,65(2):105-111.

[14] Li X S, Wang M J, Zhao X Y, et al. Finite thrust orbit interception optimization using hybrid genetic algorithm[J].ComputerSimulation, 2010, 27(12) : 20-23.(李新三,王明建,趙旭陽(yáng),等.采用混合遺傳算法的有限推力攔截軌道優(yōu)化[J].計(jì)算機(jī)仿真,2010,27(12): 20-23.)

[15] Tang Y H, Chen S L, Xu M, et al. A genetic algorithm (GA) method of orbit interception with finite thrust[J].JournalofNorthwesternPolytechnicalUniversity,2005,23(5):671-674.(湯一華,陳士櫓,徐敏,等.基于遺傳算法的有限推力軌道攔截優(yōu)化研究[J].西北工業(yè)大學(xué)學(xué)報(bào),2005,23(5):671-674.)

[16] Qi Y H, Cao X B. Design of optimal multiple-impulsive rendezvous trajectory using genetic algorithm[J].JournalofHarbinInstituteofTechnology, 2008, 40(9): 1345-1349.(齊映紅,曹喜濱.基于遺傳算法的最優(yōu)多脈沖交會(huì)軌道設(shè)計(jì)[J].哈爾濱工業(yè)大學(xué)學(xué)報(bào),2008,40(9):1345-1349.)

[17] Li Y, Li X J, Zhang D L, et al. Application of hybrid genetic algorithm in optimal Lambert orbital transfer design[J].FlightDynamics, 2013, 31(3): 269-272.(李勇,李小將,張東來(lái),等. 混合遺傳算法在最優(yōu)Lambert軌道轉(zhuǎn)移設(shè)計(jì)中的應(yīng)用[J]. 飛行力學(xué),2013,31(3):269-272.)

[18] Zhao R A.Spaceweaponorbitdesign[M]. Beijing:China Aerospace Press, 2008.(趙瑞安.空間武器軌道設(shè)計(jì)[M].北京:中國(guó)宇航出版社, 2008.)

[19] Xu P D, Huang S Y, Li Z M. A new method for determining intercepting orbit of space kinetic weapon[J].ComputerSimulation, 2009, 26(6):54-58.(徐培德,黃思勇,李志猛.天基動(dòng)能武器攔截軌道的確定方法研究[J].計(jì)算機(jī)仿真,2009,26(6):54-58.)

[20] Zhang P Y. Optimal trajectory design for space interception[D]. Harbin: Harbin Institute of Technology, 2009.(張鵬宇.空間攔截最優(yōu)軌道設(shè)計(jì)[D].哈爾濱:哈爾濱工業(yè)大學(xué),2009.)

Iterative algorithm for minimum-energy interception orbit with single pulse

MENG Shao-fei1, LIU Xin-xue1, FU Dan2, YANG Qi1

(1.PrimaryCommandCollege,TheRocketForceUniversityofEngineering,Xi’an710025,China; 2.TheFirstInstituteofChinaAerospaceScienceandTechnologyCorporation,Beijing100071,Chian)

Optimal interception orbit calculation is important content of spacecraft on-orbit interception question, and is one of the core technical problems of space-based weapons development. It’s difficult to analytically analyze and calculate because of the nonlinearity and implicitness of orbit calculation equations and interception equations. Aiming at interception trajectory optimization problem with no-constant interception point, the law of interception trajectory parameter and characteristic velocity in transfer point changing with interception time and interception true anomaly is summed up. On this basis, the calculation of the initial interception parameters from interception orbit true anomaly has been researched and rapid iterative algorithm for minimum-energy interception trajectory has been designed. Finally, practicability and rapidity of the algorithm is verified by simulation.

space interception; characteristic velocity; rapid iterative

2015-07-16;

2016-10-20;網(wǎng)絡(luò)優(yōu)先出版日期:2016-10-27。

V 448

A

10.3969/j.issn.1001-506X.2016.12.21

孟少飛(1989-),男,博士研究生,主要研究方向?yàn)轱w行動(dòng)力學(xué)與制導(dǎo)。

E-mail:627786952@qq.com

劉新學(xué)(1964-),男,教授,博士,主要研究方向?yàn)轱w行動(dòng)力學(xué)與制導(dǎo)。

E-mail:bobxaliang@yahoo.com.cn

傅 丹(1979-),女,博士,主要研究方向?yàn)轱w行器總體設(shè)計(jì)。

E-mail:1126043246@qq.com

楊 其(1983-),男,博士研究生,主要研究方向?yàn)閼T性制導(dǎo)與控制。

E-mail:epgc_yangqi@163.com

網(wǎng)絡(luò)優(yōu)先出版地址:http://www.cnki.net/kcms/detail/11.2422.TN.20161027.1600.014.html