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

孟少飛, 劉新學, 傅 丹, 楊 其

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

?

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

孟少飛1, 劉新學1, 傅 丹2, 楊 其1

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

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

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

0 引 言

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

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

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

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

(1)

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

(2)

(3)

(4)

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

1.1 特征速度關于軌道半長軸的導數(shù)

設攔截器變軌點速度增量的大小為Δv,則有

(5)

式中,v+和v-為攔截器變軌前和變軌后的速度；φ為它們的夾角。可以看出,當變軌時刻和攔截時刻都不確定時,Δv為變量v+、v-和φ的函數(shù);而當變軌時刻確定而攔截時刻不確定時,Δv僅為變量v+和φ的函數(shù)。

設a為變軌后攔截軌道的半長軸,將式(5)兩邊分別對a求導[20]有

(6)

其中

(7)

(8)

由橢圓軌道的速度公式得

(9)

對式(1)～式(9)兩邊求導可得

(10)

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

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

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

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

(11)

將式(11)兩邊對a求導可得

(12)

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

(13)

(14)

(15)

(16)

又由式(13)可得

(17)

(18)

同時,對式(1)求導可得

(19)

(20)

(21)

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

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

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

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

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

但在圖3中,攔截軌道特征速度在各區(qū)間的變化是非單調(diào)的,且攔截軌道半長軸和攔截軌道特征速度之間并非單映射的關系,無法用數(shù)值迭代的方法求解局部最優(yōu)攔截軌道。為了方便問題的解決,考慮用轉(zhuǎn)移軌道的真近點角Δθ取代轉(zhuǎn)移軌道半長軸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ī)律與變軌時刻的具體條件有關,但由于目標航天器的周期性運動,導致Δv在不同的周期內(nèi),具有相似的變化規(guī)律;在目標的一個運動周期內(nèi),當攔截軌道與原軌道方向一致時,Δv存在唯一最小值Δvmin,且以Δvmin為分界,其兩側(cè)的變化分別都是單調(diào)的。

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

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

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

則有

(24)

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

(25)

(26)

當計算精度滿足要求時,迭代停止:

|Δvn-Δvn-1|≤δ

(27)

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

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

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

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

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

(28)

(29)

(30)

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

(31)

又設

(32)

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

(33)

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

(34)

(35)

若m*=m′,則對應的ni為所求解。若m*≠m′,則對應的ni需要篩除。當n確定后,計算

(36)

(37)

設ΔθT為目標飛行過的真近點角度,則有

(38)

3 仿真算例

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

表1 目標和攔截器初始軌道根數(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

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

Δvmin=

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

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

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

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

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

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

4 結(jié) 論

文章探討了在一般天基攔截問題中,攔截終點不固定條件下的最省能量攔截軌道。首先仿真研究了變軌點速度增量隨攔截器飛過的真近點角的變化曲線并探討了其共性特點;其次在此基礎上研究設計了最省能量攔截軌道的分段迭代方法,并進行了仿真驗證,證明了該方法的有效性。此算法實現(xià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.(陳小剛,吳志林,高宇.基于能量的天基攔截軌道設計與優(yōu)化[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ù)推力最短時間轉(zhuǎn)移軌道設計[J]. 西北工業(yè)大學學報, 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)化的脈沖軌道設計[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)移的虛交點軌道攔截優(yōu)化[J].飛行力學,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]. 宇航學報,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ōu)化[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è)大學學報,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)多脈沖交會軌道設計[J].哈爾濱工業(yè)大學學報,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.(李勇,李小將,張東來,等. 混合遺傳算法在最優(yōu)Lambert軌道轉(zhuǎn)移設計中的應用[J]. 飛行力學,2013,31(3):269-272.)

[18] Zhao R A.Spaceweaponorbitdesign[M]. Beijing:China Aerospace Press, 2008.(趙瑞安.空間武器軌道設計[M].北京:中國宇航出版社, 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.(徐培德,黃思勇,李志猛.天基動能武器攔截軌道的確定方法研究[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)軌道設計[D].哈爾濱:哈爾濱工業(yè)大學,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)絡優(yōu)先出版日期:2016-10-27。

V 448

A

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

孟少飛(1989-),男,博士研究生,主要研究方向為飛行動力學與制導。

E-mail:627786952@qq.com

劉新學(1964-),男,教授,博士,主要研究方向為飛行動力學與制導。

E-mail:bobxaliang@yahoo.com.cn

傅 丹(1979-),女,博士,主要研究方向為飛行器總體設計。

E-mail:1126043246@qq.com

楊 其(1983-),男,博士研究生,主要研究方向為慣性制導與控制。

E-mail:epgc_yangqi@163.com

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