高超声速飞行器具有射程远、突防能力强和毁伤威力大等优势,已成为各国研究的热点。为准确把握高超声速飞行器的性能及特点,滑翔弹道规划技术受到广泛关注。通过求解滑翔弹道解析解,不仅可以快速规划滑翔段弹道,还可用于在线预测制导,具有重要意义。
早在20世纪50年代,滑翔弹道解析解就引起了科研人员的关注。Chapman^[1]首先给出了定升阻比下再入弹道的近似解析解;Loh^[2]则在此基础上分析了升阻比变化时的再入弹道近似解析解;Shi和Pottsepp^[3]及Naidu^[4]则采用匹配渐近展开分别求解了纵向弹道和三维滑翔弹道解析解,但在大纵程的情况下精度较低;Vinh等^[5]给出了定升阻比平衡滑翔纵向弹道解析解,具有较高的精度;Bell^[6]则进一步考虑了定升阻比和定倾侧角条件下的平衡滑翔三维弹道,但其在横向机动范围较大时误差较大。
从20世纪70年代起,随着倾侧翻转的横向机动策略的提出,高精度纵向弹道解析解成为了发展的重点。Harpold和Graves^[7]提出以D-V剖面规划再入纵向弹道,并给出了高度、高度1阶导数和射程的解析解;Yu和Chen^[8]则给出了平衡滑翔条件下定升阻比弹道的弹道倾角解析解;之后,Yu和Chen^[9]进一步给出了给定纵向升阻比和横向升阻比的三维弹道解析解;李邦杰和王明海^[10]则给出了最大升阻比平衡滑翔条件下的高度、弹道倾角、飞行时间和射程的解析解;郭兴玲和张珩^[11]及徐明亮等^[12]则给出了定弹道倾角条件下的纵向弹道解析解。
进入21世纪,随着再入飞行器横向机动能力的提高和机动突防任务的提出,滑翔段三维弹道的快速解算成为了研究的热点。Mease等^[13]推导了降阶的再入运动方程,提出了2种阻力加速度和侧向加速度混合的快速轨迹规划算法,将航天飞机的二维再入轨迹扩展到三维;Shen和Lu^[14]则利用倾侧角走廊快速规划滑翔段弹道;Rahman和Zhou等^[15-16]则通过贝塞尔曲线直接规划几何弹道,并通过调整曲线的控制参数满足约束要求。然而上述方法均需要进行弹道积分,计算量较大。
在平稳滑翔概念的基础之上^[17],本文给出了一种高精度的平稳滑翔弹道解析解。首先,将升力系数分解为横向分量、平衡滑翔纵向分量和平稳滑翔纵向分量3个部分,由此获得了相互解耦的动力学方程;然后,分别采用解析积分、正则摄动法、高斯积分法等方法对高度及射程、弹道偏角、经度及纬度进行了解析求解,并采用单步龙格-库塔积分法对速度进行求解;最后,利用上述结果,提出了一种以升力系数平稳滑翔纵向分量和横向分量为规划变量的平稳滑翔弹道快速规划方法。
1 数学模型 1.1 坐标系定义 为便于理论分析,假设地球为静止正圆球体,并引入广义赤道,如图 1所示,飞行器在广义赤道附近运动。当地地心坐标系S_e1如图 1所示,坐标原点位于地心,x_e1由地心指向飞行器;y_e1在由x_e1与速度矢量V构成的平面内,且与V同一侧;z_e1由右手定则确定。地心赤道旋转坐标系S_e坐标原点同样为地心,z_e由地心指向北极,x_e处在赤道平面内指向本初子午线,y_e由右手定则确定。两者之间的坐标转换关系为

图 1 广义赤道及坐标系 Fig. 1 Generalized equator and coordinate

图选项

(1)

式中:
l₁₁=cosωcosicosΩ-sinωsinΩ
l₁₂=cosωcosisinΩ+sinωcosΩ
l₁₃=-cosωsini
l₂₁=-sinωcosicosΩ-cosωsinΩ
l₂₂=-sinωcosisinΩ+cosωcosΩ
l₂₃=sinωsini
l₃₁=sinicosΩ
l₃₂=sinisinΩ
l₃₃=cosi
其中:i、Ω和ω分别为由地心赤道旋转坐标系向当地地心坐标系转换的轨道倾角、升交点赤经和近地幅角。
1.2 运动方程 为了便于解析计算,在当地地心坐标系下建立高超声速飞行器的动力学模型。忽略地球自转的影响,运动方程表述为

(2)

(3)

(4)

(5)

(6)

(7)

(8)

式中:h、s、θ、?、ψ和V分别为当地地心坐标系下的高度、射程、经度、纬度、弹道偏角和飞行器与地球的相对速度;γ为弹道倾角;r为地心与飞行器的距离,r=R₀+h,R₀为地球半径;m为飞行器的质量;σ为倾侧角;g为当地重力加速度,g=μ/r²,μ为地球的重力常数;L和D分别为升力及阻力,表达式为

(9)

其中:ρ为大气密度;S为气动参考面积;C_L为升力系数;C_D为阻力系数。
1.3 运动方程解耦和分段求解策略 为了便于滑翔段弹道解析求解,将升力系数分解如下:

(10)

式中:C_Y为升力系数的横向分量;C_N1为升力系数的平稳滑翔纵向分量,上述两者将作为本文平稳滑翔弹道的规划变量;C_N0为升力系数的平衡滑翔纵向分量,满足如下关系:

(11)

为了在弹道偏角解析求解时sin ?和cos ψ能够进行线性化,采用如图 2所示的分段求解策略,并且在每一段弹道的起始点均定义广义赤道坐标系(图 2中的x⁽ⁱ⁾_e1和y⁽ⁱ⁾_e1为第i段弹道起点的当地地心坐标系坐标轴),从而使得

图 2 平稳滑翔弹道的分段 Fig. 2 Segment of steady glide trajectory

图选项

(12)

同时,对于平稳滑翔的再入弹道,其弹道倾角也近似于零。进一步假定整个下滑过程中γ单调变化,则式(2)~式(8)给出的运动方程转换为如下形式:

(13)

(14)

(15)

(16)

(17)

(18)

由式(13)~式(18)可以看出,在对纵向升力系数进行分解后,纵向运动方程与其他状态变量之间完全解耦;横向运动也与速度项解耦,从而有利于解析求解的实现。
2 滑翔弹道解析解 在当地地心坐标系下,滑翔段第i段弹道的起点坐标为(0,0,h_i0);初始弹道偏角为π/2;初始速度为V_i0;初始飞行距离为0;初始弹道倾角和终端弹道倾角分别为γ_i0和γ_if,则该段弹道其他状态变量的解析解如下。
2.1 滑翔高度解析解 设C_N1为常数,对式(13)进行积分可得

(19)

式中:β_h为指数大气模型常数;ρ_i0和ρ_if分别为滑翔段第i段弹道的初始密度和终端密度,分别如下:

(20)

其中:ρ₀为海平面大气密度；h_if为第i段弹道终点高度，对式(19)整理可得

(21)

式中:K^*为纵向弹道机动常数,K^*=SC_N1ρ_i0/(2mβ_h)-cosγ_i0。在滑翔段,K^*通常大于1。
2.2 射程解析解 将式(21)代入式(14)中的路程微分方程可得

(22)

取为地心平均距离(h₀为飞行器滑翔段的初始高度),对式(22)进行积分可得

(23)

式中:s_if为第i段弹道的射程。

(24)

(25)

2.3 基于正则摄动的弹道偏角解析解 在式(21)中,K^*+1通常为10^-7数量级,而γ²/2则为10^-5数量级,因此可设K^*+cosγ≈γ²/2。在此基础上,将式(21)代入式(15)和式(16)可得

(26)

(27)

不难发现,式(26)和式(27)中,ψ与?耦合,难以直接求解。考虑到,采用正则摄动法进行估算。定义ψ和?表达式如下:

(28)

(29)

式中:ε为正则摄动常数。
由式(26)和式(28)可得

(30)

(31)

式中:A=1/ε为求解过程中的符号常数。另外,由式(29)和式(27)可得

(32)

(33)

求解式(30)即可得ψ的零阶估计为

(34)

将式(34)代入式(32),可得

(35)

式中:。对式(35)进行积分可得?的零阶估计为

(36)

将式(36)代入式(31)可得

(37)

式中:。对式(37)积分可得

(38)

将式(34)和式(38)代入式(28)可得ψ终端的估计值为

(39)

2.4 基于高斯积分公式的经度和纬度解析解 将式(21)分别代入式(16)和式(17)可得

(40)

(41)

对于横向机动能力较强的高超声速飞行器,滑翔段的ψ可能存在较大变化,若仅对cosψ和sin ψ1阶展开则会带来较大的误差。因此,将采用高斯积分公式来估算终端经度和纬度,具体如下:

(42)

(43)

式中:;;γ_j为[γ_i0,γ_if]中的第j个高斯节点;B_j为高斯求积系数。
2.5 基于单步龙格-库塔积分的滑翔速度解 将式(21)代入式(18),并取w=lnV,可得

(44)

图 3比较了V与γ的关系及w与γ的关系。图中:f_v为归一化的速度;V_c为速度归一化常数,通常取第一宇宙速度。可以看出,w与γ之间的线性化程度要远远高于V与γ之间的线性化程度。因此,式(44)可采用单步龙格-库塔积分来估算终端速度大小,具体如下:

图 3 归一化速度fv与弹道倾角γ的关系 Fig. 3 Relationship between normalized velocity fv and flight path angle γ

图选项

(45)

式(45)中的C_D可由C_L和w插值获得,而C_L的表达式如下:

(46)

最终可得积分结果如下:

(47)

式中:
w_i0=ln V_i0
w_if=ln V_if
Δγ=γ_if-γ_i0
k₁=f_K(w_i0,γ_i0)
k₂=f_K(w_i0+k₁Δγ/2,γ_i0/2+γ_if/2)
k₃=f_K(w_i0+k₂Δγ/2,γ_i0/2+γ_if/2)
k₄=f_K(w_i0+k₃Δγ,γ_if)
最后,由w_if可解得终端速度估计为

(48)

2.6 最优初始下滑角度解析解 由平稳滑翔弹道设计理论可知^[17],合理的初始弹道倾角γ₀会使整个滑翔段的控制规律更加平稳。因此,本文中初始攻角取平稳滑翔弹道倾角:

(49)

式中:V₀为滑翔段初始速度;K_N0为滑翔段初始纵向升阻比。
3 解析解精度校验 3.1 解析解与数值积分的精度对比 为了检验本文解析解的精度,采用CAV-H为飞行器模型^[18],对滑翔段弹道进行仿真计算。仿真参数设置如表 1所示。图 4~图 6分别给出了速度、纵向弹道和横向弹道的解析解与数值解的对比。可以看出,本文解析解具有较高的求解精度,其中纵向弹道的解析解与数值解几乎完全重合,而横向弹道和速度的解析解与数值解存在较小的偏差,并且随着射程的增大而增大。
表 1 仿真参数设置 Table 1 Simulation parameters setting

????????参数	数值
h0/km	60
V0/(m·s-1)	7 000
θ0/(°)	0
?0/(°)	0
ψ0/(°)	π/2
s0/km	0
hf/km	30
γ0/(°)	-0.068 4
CN1	-0.001 1
CY	0.4

表选项

图 4 速度解析解与数值解对比 Fig. 4 Comparison between velocity analytical solution and numerical solution

图选项

图 5 纵向弹道解析解与数值解对比 Fig. 5 Comparison between longitudinal trajectory analytical solution and numerical solution

图选项

图 6 横向弹道解析解与数值解对比 Fig. 6 Comparison between lateral trajectory analytical solution and numerical solution

图选项

3.2 本文解析解与Bell解析解的精度对比 将本文解析解与Bell解析解的精度进行对比。仿真的初始参数与表 1相同,终端速度设定为2 000 m/s,横向升力系数C_Y从0逐渐增加到0.7,并采用分段逐步求解的方法来预测终端状态,得到的结果如图 7~图 9所示。

图 7 弹道分为1段时的解析解精度对比 Fig. 7 Precision comparison of analytical solutions for 1 trajectory phase

图选项

图 8 弹道分为2段时的解析解精度对比 Fig. 8 Precision comparison of analytical solutions for 2 trajectory phases

图选项

图 9 弹道分为5段时的解析解精度对比 Fig. 9 Precision comparison of analytical solutions for 5 trajectory phases

图选项

由图 7可知,在采用1段弹道直接预测终端状态时,本文解析解的终端位置误差要远小于Bell解析解的终端位置误差;而本文解析解的终端速度误差则随着射程的增大逐渐增大,在射程较小时比Bell解析解精度要高,在射程较大时比Bell解析解精度要低。这是由于在给定C_N1的条件下,滑翔段弹道的C_D会有细微的变化,并且C_Y越小C_D的变化相对越明显,从而带来积分误差;而Bell解析解则是在给定纵向升阻比K_N0下获得的,其终端速度误差主要是由于忽略r的变化带来的,因此近乎为常值。由图 8和图 9可知,在采用分段逐步求解的方法后,本文解析解的精度随着步数的增加而快速提高,而Bell解析解的终端速度误差则无明显变化,终端位置误差则缓慢减小。
表 2给出了解析解精度对比。可以看出,分段数越多,本文解析解的优势越明显。当采用2段求解时,本文解析解的终端位置误差比Bell解析解的终端位置误差降低了约1个数量级,并且终端速度误差也小于Bell解析解;当采用5段求解时,本文解析解的终端位置误差则比Bell解析解的终端位置误差降低了约2个数量级,而终端速度误差则降低了约1个数量级。
表 2 解析解精度对比 Table 2 Precision comparison of analytical solutions

段数	CY	本文终端 位置误差/ km	Bell终端 位置误差/ km	本文终端 速度误差/ (m·s-1)	Bell终端 速度误差/ (m·s-1)
1	0.1	171.11	811.56	-1 940.12	-631.06
2	0.1	47.01	410.28	-642.10	-631.06
5	0.1	9.70	237.18	-44.41	-631.06
1	0.6	85.83	656.80	-90.72	-700.74
2	0.6	29.34	519.40	-41.89	-700.74
5	0.6	1.79	432.99	-3.89	-700.74

表选项






总的来说,采用本文解析解的滑翔段弹道终端状态预测精度要优于Bell解析解。
4 基于解析解的快速弹道规划 将本文解析解应用于滑翔段机动弹道规划。滑翔段初始参数与表 1相同,终端条件和过程约束见表 3,禁飞区设置见表 4。表 3中,、q_max和n_max分别为滑翔飞行器再入过程中允许的最大热流密度、最大动压和最大过载。图 10给出了应用滑翔弹道快速生成算法的滑翔段机动弹道规划流程,图中：γ_c为倾侧角反转时的弹道倾角，纵向弹道和横向弹道由式(21)、式(23)、式(39)、式(42)和式(43)确定,而终端速度约束和过程约束则由式(21)和式(48)确定。
表 3 终端条件和过程约束 Table 3 Terminal conditions and path constraints

????????参数	数值
hf/km	30
Vf/(m·s-1)	2 000
目标1经纬度/(°)	(100,0)
目标2经纬度/(°)	(115,0)
目标3经纬度/(°)	(130,0)
	500
qmax/Pa	60 000
nmax	2.5

表选项






表 4 禁飞区设置 Table 4 No-fly zone setting

禁飞区	中心坐标/(°)	半径/km
禁飞区1	(40,-6)	891
禁飞区2	(90,6)	1 002

表选项

图 10 滑翔段机动弹道规划流程 Fig. 10 Flowchart of gliding maneuvering trajectory planning

图选项

图 11~图 16给出了滑翔段机动弹道的规划结果,图 16中的EGC表示由平衡滑翔条件获得的滑翔弹道高度上界。由图 11和图 12可以看出,基于解析解规划获得的攻角α曲线和倾侧角曲线均非常光滑,并且可通过分段实现倾侧角反转。由图 13和图 14可以看出,规划弹道的终端高度和速度满足要求。由图 15可以看出,规划的滑翔段机动弹道符合禁飞区约束要求和终端经纬度要求。由图 16可以看出,规划的滑翔弹道满足最大热流密度、最大动压和最大过载等约束,并且不同任务弹道在再入走廊中几乎完全重合。

图 11 不同任务下的攻角曲线 Fig. 11 Angle of attack curves for different missions

图选项

图 12 不同任务下的倾侧角曲线 Fig. 12 Bank angle curves of different missions

图选项

图 13 不同任务下的高度曲线 Fig. 13 Altitude curves of different missions

图选项

图 14 不同任务下的速度曲线 Fig. 14 Velocity curves of different missions

图选项

图 15 不同任务下的横向弹道 Fig. 15 Lateral trajectories of different missions

图选项

图 16 不同任务下的再入走廊 Fig. 16 Reentry corridors of different missions

图选项

综上可得,利用本文的解析解能够规划出符合约束要求的滑翔弹道。在普通台式机上,在无特定初值的条件下,规划1条机动弹道的时间小于0.3 s。
5 结 论 本文给出了一种高精度的平稳滑翔弹道解析解,具有如下特点:
1) 提出了一种升力系数的分解方法,将其分解为横向分量、平衡滑翔纵向分量和平稳滑翔纵向分量3个部分,从而实现了对动力学方程中纵向运动方程、横向运动方程和速度方程的解耦。
2) 利用解耦的运动方程,首先通过解析积分获得了高度和射程的精确解析解,然后采用正则摄动法获得了较为精确的弹道偏角解析解,并在此基础上利用高斯积分公式获得了经度和纬度的解析解,最后采用单步龙格-库塔积分获得了速度的解。
3) 提出了分段逐步求解的方法来提高解析解的精度。仿真校验表明,在分段逐步求解的情况下,本文解析解的精度要大于Bell解析解的精度。
4) 采用本文解析解进行滑翔段弹道规划,不需要进行弹道积分,同时采用解析解还可先规划出符合突防要求的弹道,之后再进行速度校正,大大提高了弹道规划的效率。仿真校验表明,在普通台式机上规划1条平稳滑翔机动弹道的时间小于0.3 s。

参考文献

[1]	CHAPMAN D R.An approximate analytical method for studying entry into planetary atmospheres:NACA TN-4276[R].[S.l:s.n.],1958.
	Click to display the text
[2]	LOH W H. Some exact analytical solutions of planetary entry[J]. AIAA Journal,1963, 1(4): 836–842.
	Click to display the text
[3]	SHI Y Y, POTTSEPP L. Asymptotic expansion of a hypervelocity atmospheric entry problem[J]. AIAA Journal,1969, 7(2): 353–355.
	Click to display the text
[4]	NAIDU D S. Three-dimensional atmospheric entry problem using method of matched asymptotic expansion[J]. IEEE Transactions on Aerospace and Electronic Systems,1989, 25(5): 660–667.
	Click to display the text
[5]	VINH N X, BUSEMANN A, CULP R D. Planetary entry flight mechanics[M].Ann Arbor: University of Michigan Press, 1980: 100-138.
	Click to display the text
[6]	BELL B N.A closed-form solution to lifting reentry:AFFDL-TR-65-65[R].[S.l.:s.n.],1965.
	Click to display the text
[7]	HARPOLD J C,GRAVES C A.Shuttle entry guidance:TX.Rep.JSC-14694[R].Houston:NASA Lyndon B.Johnson Space Center,1979.
	Click to display the text
[8]	YU W,CHEN W.Guidance scheme for glide range maximization of a hypersonic vehicle:AIAA-2011-6714[R].Reston:AIAA,2011.
	Click to display the text
[9]	YU W, CHEN W. Entry guidance with real-time planning of reference based on analytical solutions[J]. Advances in Space Research,2015, 55(9): 2325–2345.
	Click to display the text
[10]	李邦杰, 王明海. 滑翔式远程导弹滑翔段弹道研究[J]. 宇航学报,2009, 30(6): 2122–2126.LI B J, WANG M H. Research on glide trajectory of long range glide missile[J]. Journal of Astronautics,2009, 30(6): 2122–2126.(in Chinese).
	Cited By in Cnki (0) | Click to display the text
[11]	郭兴玲, 张珩. 基于分段常滑翔角的长航时纵向远程滑翔飞行方程近似解[J]. 宇航学报,2005, 26(6): 712–716.GUO X L, ZHANG H. The approximate solution for long-time and long-range longitudinal gliding flight based on segment constant gliding angle[J]. Journal of Astronautics,2005, 26(6): 712–716.(in Chinese).
	Cited By in Cnki (0) | Click to display the text
[12]	徐明亮, 陈克俊, 刘鲁华, 等. 高超声速飞行器准平衡滑翔自适应制导方法[J]. 中国科学:技术科学,2012, 42(4): 378–387.XU M L, CHEN K J, LIU L H, et al. Quasi-equilibrium glide adaptive guidance for hypersonic vehicles[J]. Scientia Sinica Technologica,2012, 42(4): 378–387.(in Chinese).
	Cited By in Cnki (0) | Click to display the text
[13]	MEASE K D, CHEN D T, TEUFEL P, et al. Reduced-order entry trajectory planning for acceleration guidance[J]. Journal of Guidance,Control and Dynamics,2002, 25(2): 257–266.
	Click to display the text
[14]	SHEN Z, LU P. Onboard generation of three dimensional constrained entry trajectory[J]. Journal of Guidance,Control and Dynamic,2003, 26(1): 111–121.
	Click to display the text
[15]	ZHOU H, RAHMAN T, CHEN W. Neural network assisted inverse dynamic guidance for terminally constrained entry flight[J]. Scientific World Journal,2014, 2014(1): 171–192.
	Click to display the text
[16]	RAHMAN T,ZHOU H,CHEN W.Bezier approximation based inverse dynamic guidance for entry glide trajectory:ASCC-2013-6606111[R].Istanbul:Asian Control Conference,2013.
	Click to display the text
[17]	胡锦川, 陈万春. 高超声速飞行器平稳滑翔弹道设计方法[J]. 北京航空航天大学学报,2015, 41(8): 1464–1475.HU J C, CHEN W C. Steady glide trajectory planning method for hypersonic reentry vehicle[J]. Journal of Beijing University of Aeronautics and Astronautics,2015, 41(8): 1464–1475.(in Chinese).
	Cited By in Cnki (0) | Click to display the text
[18]	ZHANG K,CHEN W.Reentry vehicle constrained trajectory optimization:AIAA-2011-2231[R].Reston:AIAA,2011.
	Click to display the text