CN103838914A - 一种高超声速飞行器滑翔段弹道解析求解方法 - Google Patents

Abstract

本发明公开了一种高超声速飞行器滑翔段弹道解析求解方法，包括以下几个步骤：步骤1：定义当地地心坐标系；步骤2：获取在当地地心坐标系下的运动方程；步骤3：获取滑翔弹道的解析解。本发明提出了升力系数坡度率调整量的概念，使得动力学方程能够将速度项分离，获得控制量与几何弹道的直接关系；本发明提出了再入弹道分段线性化的策略，通过引入当地地心坐标系，使得再入横向弹道能够线性化。

Description

一种高超声速飞行器滑翔段弹道解析求解方法

技术领域

本发明涉及一种高超声速飞行器滑翔段弹道解析求解方法，属于航天技术、武器技术领域。 

背景技术

高超声速飞行器是指飞行马赫数大于或者等于5的飞行器，它具有飞行速度快、突防能力强、全球到达、毁伤威力大等独特的优势，已经成为当今世界各国武器研制的热点和焦点。再入弹道规划则是高超声速飞行器的一项重要的关键技术，对于典型的高超声速飞行器再入弹道，滑翔段占据航程的绝大部分，因此，滑翔段的弹道规划就显得尤为重要。 

在滑翔段的飞行过程中，飞行的速度和空域远大于常规导弹，将经历高度、马赫数的大范围变化，并且受到热流密度、动压、过载等诸多约束，使得滑翔段弹道规划成为一个高敏度问题。当前滑翔段弹道规划多采用数值优化算法来获得，如物理规划、高斯伪谱法、序列二次规划法、遗传算法+序列二次规划法、Legendre伪谱法等。然而这些方法都需要进行大量的计算，其中很大一部分是关于弹道积分的。因此，若能获得滑翔段弹道的解析解，将能够大大提高弹道规划的效率，同时，它还能够反应滑翔段飞行物理特性。因此滑翔段弹道解析解具有很高的理论意义和实际价值。 

发明内容

本发明的目的是为了解决上述问题，提出一种高超声速飞行器滑翔段弹道解析求解方法，获得高度、弹道倾角、弹道偏角、经纬度、速度、攻角和倾侧角的解析解，为快速弹道规划提供理论支持。 

一种高超声速飞行器滑翔段弹道解析求解方法，包括以下几个步骤： 

步骤1：定义当地地心坐标系； 

步骤2：获取在当地地心坐标系下的运动方程； 

步骤3：获取滑翔弹道的解析解。 

本发明的优点在于： 

（1）提出了升力系数坡度率调整量的概念，使得动力学方程能够将速度项分离，获得控制量与几何弹道的直接关系； 

（2）提出了再入弹道分段线性化的策略，通过引入当地地心坐标系，使得再入横向弹道能够线性化； 

（3）提出了以C_N1和C_Y为控制变量的再入弹道规划策略，使得横向和纵向弹道能够完全解耦； 

（4）提出了考虑纵向弹道机动的能量-射程映射关系，从而获得更精确的速度解析解。 

（5）在进行弹道规划时不需要进行弹道积分，大大加快规划运算速度； 

（6）利用解析弹道横纵解耦的特点可使得弹道规划时可以分开考虑横向弹道和纵向弹道。 

附图说明

图1是广义赤道示意图； 

图2是当地地心坐标系示意图； 

图3是地心赤道旋转坐标系与当地地心坐标系的关系图； 

图4是攻角曲线精度验证； 

图5是倾侧角曲线精度验证； 

图6是高度曲线精度验证； 

图7是横向弹道精度验证； 

图8是弹道偏角曲线精度验证； 

图9是弹道倾角曲线精度验证； 

图10是速度曲线精度验； 

图11(a)是不同初始弹道倾角下的攻角曲线； 

图11(b)是不同初始弹道倾角下的倾侧角曲线； 

图11(c)是不同初始弹道倾角下的高度曲线； 

图11(d)是不同初始弹道倾角下的横向弹道曲线； 

图11(e)是不同初始弹道倾角下的弹道倾角曲线族； 

图11(f)是不同初始弹道倾角下的H-V走廊曲线； 

图12(a)是攻角曲线族； 

图12(b)是倾侧角曲线族； 

图12(c)是高度曲线族； 

图12(d)是横向弹道曲线族； 

图12(e)是弹道偏角曲线族； 

图12(f)是弹道倾角曲线族； 

图12(g)是速度曲线族； 

图12(h)是过载曲线族； 

图12(i)是热流密度曲线族。 

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。 

本发明在拟平衡滑翔的基础上，提出了纵向升力系数余量的概念，在此基础上推导出了速度无关运动方程；利用上述方程推导获得了纵向升力系数余量与飞行距离、飞行高度和飞行弹道倾角之间的关系；在此基础上，采用分段策略，并引入当地地心坐标系，并利用它求解横向机动解析解；最后利用射程与能量之间的关系，获得了速度的解析解。 

本发明是一种高超声速飞行器滑翔段弹道解析求解方法，包括以下几个步骤： 

步骤1：定义当地地心坐标系 

为了能够更准确的求解飞行器横程，引入广义赤道的概念，如图1所示，从而使得飞行器的运动在广义赤道附近。为了描述广义赤道，引入当地地心坐标系，用符号S_e1表示，如图2所示，坐标原点位于地心；x_e1由地心指向飞行器；y_e1在由x_e1与V的平面内，且与V同一侧；z_e1由右手定则确定(地心赤道旋转坐标系坐标原点同样为地心，z_e由地心北极，x_e处在赤道平面内指向本初子午线，y_e由右手定则确定)。 

进一步，地心赤道旋转坐标系与当地地心坐标系的关系，如图3所示。其中：i₁为轨道倾角（所谓轨道，即指广义赤道），取值范围为0°至180°；Ω₁为升交点赤经，取值范围为0°至360°；ω₁为近地点幅角，取值范围为0°至360°。x1、y1、y2均为过程坐标 轴，用于描述坐标转换的过程。 

（1）求解转换矩阵 

设飞行器的当前经纬度为[θ,φ]，则它在地心赤道旋转坐标系中的单位位置矢量为， 

$\stackrel{\→}{r}={\left[\begin{array}{ccc}\cos \mathrm{\φ}\cos \mathrm{\θ}& \cos \mathrm{\φ}\sin \mathrm{\θ}& \sin \mathrm{\φ}\end{array}\right]}^T---\left(1\right)$

其中：θ为在地心赤道旋转坐标系下的当地经度；φ为地心赤道旋转坐标系下的纬度；

为飞行器当前位置矢量。 

为了在地心赤道旋转坐标系中将单位速度矢量表示出，需要进行如下坐标转换， 

其中：S_e表示地心赤道旋转坐标系；S_u为当地铅垂坐标系；S_k为航迹坐标系。 

ψ为地心赤道旋转坐标系下的行向角；γ为地心赤道旋转坐标系下的弹道倾角，最终，可以将在[θ,φ]处且弹道倾角为ψ的单位速度矢量表述为， 

$\stackrel{\→}{v}=\left[\begin{array}{c}-\sin \mathrm{\ψ}\sin \mathrm{\θ}-\cos \mathrm{\ψ}\sin \mathrm{\φ}\cos \mathrm{\θ}\\ \sin \mathrm{\ψ}\cos \mathrm{\θ}-\cos \mathrm{\ψ}\sin \mathrm{\φ}\sin \mathrm{\θ}\\ \cos \mathrm{\ψ}\cos \mathrm{\φ}\end{array}\right]---\left(3\right)$

其中：

为在地心赤道旋转坐标系下的飞行器速度矢量。 

利用式(1)和式(3)，可得广义赤道平面的单位法向量为： 

${\stackrel{\→}{z}}_{e1}=\frac{\stackrel{\→}{r}\×\stackrel{\→}{v}}{|\stackrel{\→}{r}\×\stackrel{\→}{v}|}=\left[\begin{array}{c}\sin \mathrm{\θ}\cos \mathrm{\ψ}-\sin \mathrm{\φ}\sin \mathrm{\ψ}\cos \mathrm{\θ}\\ -\cos \mathrm{\θ}\cos \mathrm{\ψ}-\sin \mathrm{\φ}\sin \mathrm{\ψ}\sin \mathrm{\θ}\\ \cos \mathrm{\φ}\sin \mathrm{\ψ}\end{array}\right]---\left(4\right)$

其中：

为广义赤道平面单位法向量。而地心赤道旋转坐标系z轴的单位向量为，  ${\stackrel{\→}{z}}_e={\left[\begin{array}{ccc}0& 0& 1\end{array}\right]}^T,$ 由轨道倾角i₁的定义可知 $\cos i_1={\stackrel{\→}{z}}_{e1}^T{\stackrel{\→}{z}}_e,$ 结合式(4)可得i₁的取值为， 

i₁＝arccos(cosφsinψ)          (5) 

利用

和

可进一步获得当地地心坐标系的x轴单位向量

${\stackrel{\→}{x}}_1=\frac{{\stackrel{\→}{z}}_e\×{\stackrel{\→}{z}}_{e1}}{|{\stackrel{\→}{z}}_e\×{\stackrel{\→}{z}}_{e1}|}=\frac{1}{\sqrt{{\cos }^2\mathrm{\ψ}+{\sin }^2\mathrm{\φ}{\sin }^2\mathrm{\ψ}}}\left[\begin{array}{c}\cos \mathrm{\θ}\cos \mathrm{\ψ}+\sin \mathrm{\φ}\sin \mathrm{\ψ}\sin \mathrm{\θ}\\ \sin \mathrm{\θ}\cos \mathrm{\ψ}-\sin \mathrm{\φ}\sin \mathrm{\ψ}\cos \mathrm{\θ}\\ 0\end{array}\right]---\left(6\right)$

上式在ψ＝π/2,φ＝0时存在歧义，此时，

由于 ${\stackrel{\→}{x}}_e=\left[\mathrm{1,0,0}\right],$ 并且升交点赤经Ω₁满足 ${\cos \mathrm{\Ω}}_1={\stackrel{\→}{x}}_1^T{\stackrel{\→}{x}}_e,$ 利用 ${\stackrel{\→}{x}}_{e1}$ 和 ${\stackrel{\→}{x}}_e$ 可获得Ω₁如下所示： 

${\mathrm{\&Omega;}}_1=\left\{\begin{array}{cc}\arccos \left(\frac{\cos \mathrm{\&theta;}\cos \mathrm{\&psi;}+\sin \mathrm{\&phi;}\sin \mathrm{\&psi;}\sin \mathrm{\&theta;}}{\sqrt{{\cos }^2\mathrm{\&psi;}+{\sin }^2\mathrm{\&phi;}{\sin }^2\mathrm{\&psi;}}}\right)& \frac{\sin \mathrm{\&theta;}\cos \mathrm{\&psi;}-\sin \mathrm{\&phi;}\sin \mathrm{\&psi;}\cos \mathrm{\&theta;}}{\sqrt{{\cos }^2\mathrm{\&psi;}{+\sin }^2\mathrm{\&phi;}{\sin }^2\mathrm{\&psi;}}}\&GreaterEqual;0\\ 2\mathrm{\&pi;}-\arccos \left(\frac{\cos \mathrm{\&theta;}\cos \mathrm{\&psi;}+\sin \mathrm{\&phi;}\sin \mathrm{\&psi;}\sin \mathrm{\&theta;}}{\sqrt{{\cos }^2\mathrm{\&psi;}+{\sin }^2\mathrm{\&phi;}{\sin }^2\mathrm{\&psi;}}}\right)& \frac{\sin \mathrm{\&theta;}\cos \mathrm{\&psi;}-\sin \mathrm{\&phi;}\sin \mathrm{\&psi;}\cos \mathrm{\&theta;}}{\sqrt{{\cos }^2\mathrm{\&psi;}+{\sin }^2\mathrm{\&phi;}{\sin }^2\mathrm{\&psi;}}}<0\end{array}\right.---\left(7\right)$

由于， ${\stackrel{\&RightArrow;}{x}}_{e1}=\stackrel{\&RightArrow;}{r},$ 并且 ${\cos \mathrm{\&omega;}}_1={\stackrel{\&RightArrow;}{x}}_1^T{\stackrel{\&RightArrow;}{x}}_{e1},$ 利用

和

可得， 

$\mathrm{judge}={\stackrel{\&RightArrow;}{z}}_{e1}^T({\stackrel{\&RightArrow;}{x}}_1\&times;{\stackrel{\&RightArrow;}{x}}_{e1})=\frac{\sin \mathrm{\&phi;}}{{\cos }^2\mathrm{\&psi;}+{\sin }^2\mathrm{\&phi;}{\sin }^2\mathrm{\&psi;}}---\left(8\right)$

${\mathrm{\&omega;}}_1=\left\{\begin{array}{cc}\arccos \left(\frac{\cos \mathrm{\&phi;}\cos \mathrm{\&psi;}}{\sqrt{{\cos }^2\mathrm{\&psi;}+{\sin }^2\mathrm{\&phi;}{\sin }^2\mathrm{\&psi;}}}\right)& \mathrm{judge}>0\\ 2\mathrm{\&pi;}-\arccos \left(\frac{\cos \mathrm{\&phi;}\cos \mathrm{\&psi;}}{\sqrt{{\cos }^2\mathrm{\&psi;}+{\sin }^2\mathrm{\&phi;}{\sin }^2\mathrm{\&psi;}}}\right)& \mathrm{judge}\&le;0\end{array}\right.---\left(9\right)$

其中：ω₁为近地点幅角，judge为判断近地点幅角正负号的特征量。 

最终，获得S_e1与S_e之间的转换关系如下所示。 

$\begin{array}{c}{L}_{e_{1}e}=L_z\left({\mathrm{\&omega;}}_1\right)L_x\left(i_1\right)L_z\left({\mathrm{\&Omega;}}_1\right)\\ =\left[\begin{array}{ccc}\left(\begin{array}{c}\cos {\mathrm{\&omega;}}_1\cos i_1\cos {\mathrm{\&Omega;}}_1\\ -\sin {\mathrm{\&omega;}}_1\sin {\mathrm{\&Omega;}}_1\end{array}\right)& \left(\begin{array}{c}\cos {\mathrm{\&omega;}}_1\cos i_1\sin {\mathrm{\&Omega;}}_1\\ +\sin {\mathrm{\&omega;}}_1\cos {\mathrm{\&Omega;}}_1\end{array}\right)& -\cos {\mathrm{\&omega;}}_1\sin i_1\\ \left(\begin{array}{c}-\sin {\mathrm{\&omega;}}_1\cos i_1\cos {\mathrm{\&Omega;}}_1\\ -\cos {\mathrm{\&omega;}}_1\sin {\mathrm{\&Omega;}}_1\end{array}\right)& \left(\begin{array}{c}-\sin {\mathrm{\&omega;}}_1\cos i_1\sin {\mathrm{\&Omega;}}_1\\ +\cos {\mathrm{\&omega;}}_1\cos {\mathrm{\&Omega;}}_1\end{array}\right)& \sin {\mathrm{\&omega;}}_1\sin i_1\\ \sin i_1\cos {\mathrm{\&Omega;}}_1& \sin i_1\sin {\mathrm{\&Omega;}}_1& \cos i_1\end{array}\right]\end{array}----\left(11\right)$

其中：

表示当地地心坐标系与地心赤道旋转坐标系之间的转换矩阵；L_z(ω₁)表示绕z轴旋转ω₁对应的转换矩阵；L_x(i₁)表示绕x轴旋转i₁对应的转换矩阵；L_z(Ω₁)表示绕z轴旋转Ω₁对应的转换矩阵。由于线性化的动力学模型存在于当地地心坐标系中，而最终获得的弹道结果需要在地心赤道旋转坐标系中表示，因此，需要

来完成两者之间的转换。 

（2）坐标转换 

设在当地地心坐标系中存在一点[θ₁,φ₁]，并且该点的弹道偏角为ψ₁，假设该点在地心赤 道旋转坐标系中的坐标为[θ,φ]，弹道偏角为ψ，则它们满足， 

$\left[\begin{array}{c}\cos \mathrm{\&phi;}\cos \mathrm{\&theta;}\\ \cos \mathrm{\&phi;}\sin \mathrm{\&theta;}\\ \sin \mathrm{\&phi;}\end{array}\right]={\left(L_{e_{1}e}\right)}^T\left[\begin{array}{c}\cos {\mathrm{\&phi;}}_1\cos {\mathrm{\&theta;}}_1\\ \cos {\mathrm{\&phi;}}_1\sin {\mathrm{\&theta;}}_1\\ \sin {\mathrm{\&phi;}}_1\end{array}\right]---\left(12\right)$

从而，获得在地心赤道旋转坐标系中的位置为， 

$\begin{array}{c}\mathrm{\&phi;}=\arcsin (l_{11}\cos {\mathrm{\&phi;}}_1\cos {\mathrm{\&theta;}}_1+l_{21}\cos {\mathrm{\&phi;}}_1\sin {\mathrm{\&theta;}}_1+l_{31}\sin {\mathrm{\&phi;}}_1)\\ \begin{array}{}\end{array},\mathrm{\&theta;}=\left\{\begin{array}{cc}\arcsin \left(\frac{l_{12}\cos {\mathrm{\&phi;}}_1\cos {\mathrm{\&theta;}}_1+l_{22}\cos {\mathrm{\&phi;}}_1\sin {\mathrm{\&theta;}}_1+l_{32}\sin {\mathrm{\&phi;}}_1}{\cos \mathrm{\&phi;}}\right)& k_{\mathrm{\&theta;}}>0\\ \arcsin \left(\frac{l_{12}\cos {\mathrm{\&phi;}}_1\cos {\mathrm{\&theta;}}_1+l_{22}\cos {\mathrm{\&phi;}}_1\sin {\mathrm{\&theta;}}_1+l_{32}\sin {\mathrm{\&phi;}}_1}{\cos \mathrm{\&phi;}}\right)+\mathrm{\&pi;}& k_{\mathrm{\&theta;}}\&le;0\end{array}\right.\end{array}---\left(13\right)$

其中： 

l₁₁＝cosω₁cosi₁cosΩ₁-sinω₁sinΩ₁

l₂₁＝-sinω₁cosi₁cosΩ₁-cosω₁sinΩ₁

l₃₁＝sini₁cosΩ₁

l₁₂＝cosω₁cosi₁sinΩ₁+sinω₁cosΩ₁

l₂₂＝-sinω₁cosi₁sinΩ₁+cosω₁cosΩ₁

l₃₂＝sini₁sinΩ₁

k_θ＝l₁₁cosφ₁cosθ₁+l₂₁cosφ₁sinθ₁+l₃₁sinφ₁

单位速度矢量定义如下， 

$\left[\begin{array}{c}{v}_x\\ v_y\\ v_z\end{array}\right]={\left(L_{e_{1}e}\right)}^T\left[[\begin{array}{c}-\sin {\mathrm{\&psi;}}_1\sin {\mathrm{\&theta;}}_1-\cos {\mathrm{\&psi;}}_1\sin {\mathrm{\&phi;}}_1\cos {\mathrm{\&theta;}}_1\\ \sin {\mathrm{\&psi;}}_1\cos {\mathrm{\&theta;}}_1-\cos {\mathrm{\&psi;}}_1\sin {\mathrm{\&phi;}}_1\sin {\mathrm{\&theta;}}_1\\ \cos {\mathrm{\&psi;}}_1\cos {\mathrm{\&phi;}}_1\end{array}\right]---\left(14\right)$

另外，在地心赤道旋转坐标系中可得， 

$\left[\begin{array}{c}-\sin \mathrm{\&psi;}\sin \mathrm{\&theta;}-\cos \mathrm{\&psi;}\sin \mathrm{\&phi;}\cos \mathrm{\&theta;}\\ \sin \mathrm{\&psi;}\cos \mathrm{\&theta;}-\cos \mathrm{\&psi;}\sin \mathrm{\&phi;}\sin \mathrm{\&theta;}\\ \cos \mathrm{\&psi;}\cos \mathrm{\&phi;}\end{array}\right]=\left[\begin{array}{c}{v}_x\\ v_y\\ v_z\end{array}\right]$

由式(14)可以解得， 

$\begin{array}{c}\cos \mathrm{\&psi;}=\frac{v_z}{\cos \mathrm{\&phi;}}\\ \sin \mathrm{\&psi;}=\left\{\begin{array}{cc}\frac{v_x+\cos \mathrm{\&psi;}\sin \mathrm{\&phi;}\cos \mathrm{\&theta;}}{-\sin \mathrm{\&theta;}}& 0\&Element;\{0,\mathrm{\&pi;}\}\\ \frac{v_y+\cos \mathrm{\&psi;}\sin \mathrm{\&phi;}\sin \mathrm{\&theta;}}{\cos \mathrm{\&theta;}}& 0\&Element;(0,\mathrm{\&pi;})\&cup;(\mathrm{\&pi;},2\mathrm{\&pi;})\end{array}\right.\end{array}---\left(15\right)$

从而可得， 

$\mathrm{\&psi;}=\left\{\begin{array}{cc}\arccos \left(\frac{v_z}{\cos \mathrm{\&phi;}}\right)& \sin \mathrm{\&psi;}\&GreaterEqual;0\\ 2\mathrm{\&pi;}-\arccos \left(\frac{v_z}{\cos \mathrm{\&phi;}}\right)& \sin \mathrm{\&psi;}<0\end{array}\right.---\left(16\right)$

其中：v_y为单位速度矢量y轴分量；v_z为单位速度矢量z轴分量。 

步骤2：获取在当地地心坐标系下的运动方程 

无动力滑翔飞行器单位质量机械能表述为e＝μ/r-V²/2，忽略地球自转，在当地地心坐标系中三自由度运动模型可表述为， 

$\begin{array}{c}\stackrel{\&CenterDot;}{h}=V\sin \mathrm{\&gamma;}\\ {\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_1=\frac{V\cos \mathrm{\&gamma;}\sin {\mathrm{\&psi;}}_1}{(h+R_0)\cos {\mathrm{\&phi;}}_1}\\ {\stackrel{\&CenterDot;}{\mathrm{\&phi;}}}_1=\frac{V\cos \mathrm{\&gamma;}\cos {\mathrm{\&psi;}}_1}{h+R_0}\\ \stackrel{\&CenterDot;}{e}=\frac{\mathrm{DV}}{m}\\ \stackrel{\&CenterDot;}{\mathrm{\&gamma;}}=\frac{1}{V}[\frac{L\cos \mathrm{\&sigma;}}{m}+(\frac{V^2}{h+R_0}-g)\cos \mathrm{\&gamma;}]\\ {\stackrel{\&CenterDot;}{\mathrm{\&psi;}}}_1=\frac{1}{V}[\frac{L\sin \mathrm{\&sigma;}}{m\cos \mathrm{\&gamma;}}+\frac{V^2}{h+R_0}\cos \mathrm{\&gamma;}\sin {\mathrm{\&psi;}}_1\tan {\mathrm{\&phi;}}_1]\\ \stackrel{\&CenterDot;}{s}=\frac{V\cos \mathrm{\&gamma;}}{h+R_0}\end{array}---\left(17\right)$

其中：R₀＝6378km，θ₁和φ₁为当地地心坐标系下的经度和纬度；V为飞行器与地球的相对速度；ψ₁为当地地心坐标系下的弹道偏角；γ为地心赤道旋转坐标系下的弹道倾角(与当地地心左边系下的弹道倾角γ₁相等)；h为地心赤道旋转坐标系下的高度（与当地地心坐标系下的高度相等）；s为滑翔段的飞行路程；

和

分别为高度、当地地心坐标系下的经度、当地地心坐标系下的纬度、负比能量、弹道倾角和当地地心坐标系下的弹道偏角及射程对时间的一阶导数。D和L分别为以过载形式表示的升力和阻力，即D＝0.5ρV²SC_D，L＝0.5ρV²SC_L，其中ρ为大气密度，S为气动参考面积，m为飞行器的质量；C_L和C_D分别为升力系数和阻力系数；g为当地重力加速度；σ为速度滚转角。 

由于φ₁≈0、γ≈0、ψ₁≈π/2，并假设0.5ρV²SC_N0＝-m(V²/r-g)cosγ且 

C_Lcosσ＝C_N0+C_N1，其中C_N0为纵向平衡滑翔升力系数、C_N1为升力系数坡度率调整量，式(17)可进一步化为： 

$\begin{array}{c}\begin{array}{cc}\stackrel{\&CenterDot;}{h}=V\sin & \stackrel{\&CenterDot;}{e}=\frac{\mathrm{DV}}{m}\\ {\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_1=\frac{V}{h+R_0}& \stackrel{\&CenterDot;}{\mathrm{\&gamma;}}=\frac{\mathrm{\&rho;}{\mathrm{VSC}}_{N1}}{2m}\\ {\stackrel{\&CenterDot;}{\mathrm{\&phi;}}}_1=\frac{V(\mathrm{\&pi;}/2-{\mathrm{\&psi;}}_1)}{h+R_0}& {\stackrel{\&CenterDot;}{\mathrm{\&psi;}}}_1=\frac{\mathrm{\&rho;}{\mathrm{VSC}}_Y}{2m}+\frac{V}{h+R_0}{\mathrm{\&phi;}}_1\end{array}\\ \stackrel{\&CenterDot;}{s}=\frac{V\cos \mathrm{\&gamma;}}{h+R_0}\end{array}---\left(18\right)$

其中：C_N1为升力系数坡度率调整量，C_Y为升力系数横向调整量，满足C_Y＝C_Lsinσ。 

取ψ₂＝π/2-ψ₁，并对上式进行如下的变换， 

$\frac{\mathrm{dh}}{\sin \mathrm{\&gamma;}}=\mathrm{Vdt}---\left(19\right)$

(h+R₀)dθ₁＝Vdt            (20) 

$\frac{(h+R_0)d{\mathrm{\&phi;}}_1}{{\mathrm{\&psi;}}_2}=\mathrm{Vdt}---\left(21\right)$

$\frac{2\mathrm{md\&gamma;}}{\mathrm{\&rho;}{\mathrm{SC}}_{N1}}=\mathrm{Vdt}---\left(22\right)$

$\frac{d{\mathrm{\&psi;}}_2}{-\frac{\mathrm{\&rho;}{\mathrm{SC}}_Y}{2m}-\frac{1}{h+R_0}{\mathrm{\&phi;}}_1}=\mathrm{Vdt}---\left(23\right)$

$\frac{\mathrm{mde}}{D}=\mathrm{Vdt}---\left(24\right)$

(h+R₀)ds＝Vdt            (25) 

其中：dh、dθ₁、dφ₁、de、dγ、dψ₁和ds分别为高度、当地地心坐标系下的经度、当地地心坐标系下的纬度、负比能量、弹道倾角和当地地心坐标系下的弹道偏角及射程的微分，dt为时间的微分。 

步骤3：滑翔弹道解析解 

设起滑点坐标为(θ₀,φ₀)，起滑高度为h₀，起滑弹道倾角为γ₀，终端高度为h_f，纵向升力系数余量为C_N1，横向升力系数为C_Y。利用式(11)可获得该点处当地地心坐标系与地心赤道 旋转坐标系之间的转换关系，并且该点在当地地心坐标系中满足θ₁₀＝0、φ₁₀＝0以及ψ₁₀＝π/2。则在当地地心坐标系下滑翔段弹道的解析解如下。 

（1）滑翔弹道倾角解析解 

由式(19)和式(22)可得， 

$\frac{\mathrm{dh}}{\sin \mathrm{\&gamma;}}=\frac{2\mathrm{md\&gamma;}}{\mathrm{\&rho;}{\mathrm{SC}}_{N1}}---\left(26\right)$

对式(26)进行积分可得， 

${\mathrm{\&gamma;}}_f-\cos {\mathrm{\&gamma;}}_0=\frac{{\mathrm{SC}}_{N1}}{2m{\mathrm{\&beta;}}_h}({\mathrm{\&rho;}}_f-{\mathrm{\&rho;}}_0)---\left(27\right)$

其中：ρ_f为预测终点大气密度；γ_f为预测终点弹道倾角。β_h为指数大气模型常数。对式(27)进一步处理可得， 

$\frac{{\mathrm{SC}}_{N1}}{2m{\mathrm{\&beta;}}_h}{\mathrm{\&rho;}}_f-\cos {\mathrm{\&gamma;}}_f=\frac{{\mathrm{SC}}_{N1}}{2m{\mathrm{\&beta;}}_h}{\mathrm{\&rho;}}_0-\cos {\mathrm{\&gamma;}}_0---\left(28\right)$

从而可以得到一个对滑翔段弹道影响非常大的常数， 

$K^{*}=\frac{{\mathrm{SC}}_{N1}{\mathrm{\&rho;}}_0}{2m{\mathrm{\&beta;}}_h}-\cos {\mathrm{\&gamma;}}_0---\left(29\right)$

上述的K^*即为决定滑翔段纵向弹道弯曲程度的常数。 

（2）飞行路程解析解 

由式(20)和式(22)可得： 

$(h+R_0)\mathrm{ds}=\frac{2\mathrm{md\&gamma;}}{\mathrm{\&rho;}{\mathrm{SC}}_{N1}}---\left(30\right)$

将式(28)和式(29)带入式(30)可得： 

${\mathrm{\&beta;}}_h(h+R_0)\mathrm{ds}=\frac{1}{\cos \mathrm{\&gamma;}+K^{*}}\mathrm{d\&gamma;}---\left(31\right)$

设

为地心平均距离。对式(31)进行积分可得： 

${\mathrm{\&beta;}}_h\stackrel{\&OverBar;}{R}{s}_f=\left\{\begin{array}{cc}\frac{1}{1+K^{*}}\sqrt{\frac{1+K^{*}}{1-K^{*}}}\left(\begin{array}{c}\ln \left|\frac{\tan \frac{{\mathrm{\&gamma;}}_f}{2}+\sqrt{\frac{1+K^{*}}{1-K^{*}}}}{\tan \frac{{\mathrm{\&gamma;}}_f}{2}-\sqrt{\frac{1+K^{*}}{1-K^{*}}}}\right|\\ -\ln |\frac{\tan \frac{{\mathrm{\&gamma;}}_0}{2}+\sqrt{\frac{1+K^{*}}{1-K^{*}}}}{\tan \frac{{\mathrm{\&gamma;}}_0}{2}-\sqrt{\frac{1+K^{*}}{1-K^{*}}}}\end{array}\right)& \left|K^{*}\right|<1\\ \frac{2}{(1+K^{*})}\sqrt{\frac{1+K^{*}}{K^{*}-1}}\left(\begin{array}{c}a\tan \left(\sqrt{\frac{K^{*}-1}{K^{*}+1}}\tan \frac{{\mathrm{\&gamma;}}_f}{2}\right)\\ -a\tan \left(\sqrt{\frac{K^{*}-1}{K^{*}+1}}\tan \frac{{\mathrm{\&gamma;}}_0}{2}\right)\end{array}\right)& \left|K^{*}\right|>1\\ \frac{2}{K^{*}-1}(\cot \frac{{\mathrm{\&gamma;}}_f}{2}-\cot \frac{{\mathrm{\&gamma;}}_0}{2})& \left|K^{*}\right|=1\end{array}\right.---\left(32\right)$

其中：s_f为预测终点航程（初始航程为0）。 

（3）滑翔经度解析解 

由式(20)和式(22)可得， 

$(h+R_0)d{\mathrm{\&theta;}}_1=\frac{2\mathrm{md\&gamma;}}{\mathrm{\&rho;}{\mathrm{SC}}_{N1}}---\left(33\right)$

式(33)与式(30)比较可得， 

dθ₁＝ds        (34) 

从而可得， 

θ_1f＝s_f            (35) 

其中：θ_1f为当地地心坐标系下的预测终点经度。 

（4）滑翔纬度与弹道偏角解析解 

由式(21)、式(22)和式(23)可得， 

$\begin{array}{c}\frac{d{\mathrm{\&phi;}}_1}{\mathrm{d\&gamma;}}=\frac{2m}{(h+R_0)\mathrm{\&rho;}{\mathrm{SC}}_{N1}}{\mathrm{\&psi;}}_2\\ \frac{d{\mathrm{\&psi;}}_2}{\mathrm{d\&gamma;}}=-\frac{C_Y}{C_{N1}}{-}_{}\frac{2m}{(h+R_0){\mathrm{\&rho;SC}}_{N1}}{\mathrm{\&phi;}}_1\end{array}---\left(36\right)$

将式(28)和式(29)带入式(36)，可得， 

$\begin{array}{c}\frac{d{\mathrm{\&phi;}}_1}{\mathrm{d\&gamma;}}=\frac{{\mathrm{\&psi;}}_2}{\stackrel{\&OverBar;}{R}{\mathrm{\&beta;}}_h\left(K^{*}+\cos \mathrm{\&gamma;}\right)}\\ \frac{\mathrm{d\&psi;}{2}_{}}{\mathrm{d\&gamma;}}=-\frac{C_Y}{C_{N1}}-\frac{{\mathrm{\&phi;}}_1}{\stackrel{\&OverBar;}{R}{\mathrm{\&beta;}}_h(K^{*}+\cos \mathrm{\&gamma;})}\end{array}---\left(37\right)$

由于φ₁₀＝0，并且在整段弹道中φ₁为小量，故

对ψ₂的影响不大，在这里取， 

$\begin{array}{c}{K}_1=\frac{1}{{\mathrm{\&gamma;}}_f-{\mathrm{\&gamma;}}_0}{\&Integral;}_{{\mathrm{\&gamma;}}_0}^{{\mathrm{\&gamma;}}_f}\frac{1}{\stackrel{\&OverBar;}{R}{\mathrm{\&beta;}}_h(K^{*}+\cos \mathrm{\&gamma;})}\mathrm{d\&gamma;}\\ =\frac{{\mathrm{\&theta;}}_{1f}}{{\mathrm{\&gamma;}}_f-{\mathrm{\&gamma;}}_0}\end{array}---\left(38\right)$

则，式(37)可化为， 

$\begin{array}{c}\frac{d{\mathrm{\&phi;}}_1}{\mathrm{d\&gamma;}}=K_1{\mathrm{\&psi;}}_2\\ \frac{d{\mathrm{\&psi;}}_2}{\mathrm{d\&gamma;}}=-\frac{C_Y}{C_{N1}}-K_1{\mathrm{\&phi;}}_1\end{array}---\left(39\right)$

由于φ₁₀＝0，ψ₂₀＝0，对上式进行拉式变换可得， 

$\begin{array}{c}s{\mathrm{\&phi;}}_1\left(s\right)=K_1{\mathrm{\&psi;}}_2\left(s\right)\\ s{\mathrm{\&psi;}}_2\left(s\right)=-\frac{C_Y}{C_{N1}}\frac{1}{s}-K_1{\mathrm{\&phi;}}_1\left(s\right)\end{array}---\left(40\right)$

由上式可得， 

$\begin{array}{c}{\mathrm{\&psi;}}_2\left(s\right)=\frac{-\frac{C_Y}{C_{N1}}}{(s^2+K_1^2)}\\ {\mathrm{\&phi;}}_1\left(s\right)=\frac{-\frac{C_Y}{C_{N1}}{K}_1}{s(s^2+K_1^2)}\end{array}---\left(41\right)$

对式(41)进行拉式反变换可得， 

$\begin{array}{c}{\mathrm{\&psi;}}_{2f}=-\frac{C_Y}{C_{N1}}\frac{1}{K_1}\sin {\mathrm{\&theta;}}_{1f}\\ {\mathrm{\&phi;}}_{1f}=-\frac{C_Y}{C_{N1}}\frac{1}{K_1}(1-\cos {\mathrm{\&theta;}}_{1f})\end{array}---\left(42\right)$

其中：φ_1f为当地地心坐标系下的预测终点纬度。由上式还可得当地地心坐标系下的预测终点弹道偏角ψ_1f为ψ_1f＝π/2-ψ_2f。利用式(42)和式(12)给出的坐标转换关系，就可以获得在 地心赤道旋转坐标系下预测终点经度θ_f预测终点纬度φ_f和预测终点弹道偏角ψ_f。 

（5）滑翔速度解析解 

假定在滑翔过程中，保持定常纵向升阻比K_N（K_N＝C_Lcosσ/C_D），则滑翔阻力可写为， 

$D=\frac{L\cos \mathrm{\&sigma;}}{K_N}---\left(43\right)$

由于Lcosσ＝0.5ρV²S(C_N0+C_N1)，而C_N1为一小量，忽略后可得， 

$D=\frac{m[g-V^2/(R_0+h)]\cos \mathrm{\&gamma;}}{K_N}+\frac{0.5\mathrm{\&rho;}{V}^2{\mathrm{SC}}_{N1}}{K_N}---\left(44\right)$

其中：D为阻力；m为飞行器质量；V为飞行器相对地球的速度；C_N1为升力系数坡度率调整量；g为当地重力加速度，还可以写成如下形式， 

$g=\frac{\mathrm{\&mu;}}{{(R_0+h)}^2}$

其中：μ为由引力常数和地球质量决定的常数。 

式(44)中，阻力可分为两个部分，即保持平衡滑翔的升力产生的阻力以及用于纵向转弯的升力产生的阻力，其中第一项占主要部分，第二项仅为小的修正。忽略第二部分的影响，第一部分阻力满足， 

$\frac{\mathrm{de}}{d{\mathrm{\&theta;}}_1}=\frac{-\frac{\mathrm{\&mu;}}{R_0+h}+2e}{K_N}---\left(45\right)$

由于h相对R₀来说是小量，这里用它的算术平均值代替，

则上式积分可得， 

$\ln \left(\frac{2e-\mathrm{\&mu;}/\stackrel{\&OverBar;}{R}}{2e_0-\mathrm{\&mu;}/\stackrel{\&OverBar;}{R}}\right)=\frac{2{\mathrm{\&theta;}}_{1f}}{K_N}---\left(46\right)$

其中：e₀为初始负比能量。从而可得， 

$2e_f=\mathrm{\&mu;}/\stackrel{\&OverBar;}{R}+(2e_0-\mathrm{\&mu;}/\stackrel{\&OverBar;}{R})\exp \left(\frac{2{\mathrm{\&theta;}}_{1f}}{K}\right)---\left(47\right)$

其中：e_f为预测终点负比能量。对于另一部分用于转弯的升力，其消耗的机械能为， 

$\mathrm{\&Delta;e}={\&Integral;}_{s_0}^s\frac{V\stackrel{\&CenterDot;}{\mathrm{\&gamma;}}}{K_N}\mathrm{rds}={\&Integral;}_t^{t_f}\frac{V^2\stackrel{\&CenterDot;}{\mathrm{\&gamma;}}}{K_N}\mathrm{dt}---\left(48\right)$

其中：为弹道倾角的导数。 

取 ${\stackrel{\&OverBar;}{V}}^2=(V_0^2+V_f^2+V_0{V}_f)/3,$ 则可得， 

$\mathrm{\&Delta;e}=\frac{{\stackrel{\&OverBar;}{V}}^2}{K_N}({\mathrm{\&gamma;}}_f-{\mathrm{\&gamma;}}_0)---\left(49\right)$

其中：Δe为负比能量修正量，综上可得，终端速度预测为， 

$V_f=\sqrt{\mathrm{\&mu;}/\stackrel{\&OverBar;}{R}-(2e_0-\mathrm{\&mu;}/\stackrel{\&OverBar;}{R})\exp (2{\mathrm{\&theta;}}_{1f}/K_N)-\frac{2{\stackrel{\&OverBar;}{V}}^2}{K_N}({\mathrm{\&gamma;}}_f-{\mathrm{\&gamma;}}_0)}---\left(50\right)$

其中：V_f为预测终点速度；

为距离地心的平均距离；K_N为升阻比。需要注意的是，在式(50)的计算中利用了平均速度，因此需要通过迭代来不断修正

同时，还可以利用迭代修正对K_N的估计（在实际弹道中K_N不为常值，为了估算的精确一般采用K_N的平均值，为了简单起见，可令

其中K_N0为起始点处的纵向升阻比，）

采取迭代修正的方式获得。 

（6）滑翔攻角和倾侧角解析解 

式(28)、式(35)、式(42)、式(50)分别给出了预测终点弹道倾角γ_f、经度θ_f、纬度φ_f、弹道偏角ψ_f、速度V_f与高度h_f之间的关系。利用平衡滑翔关系可知， 

$C_{N0}=-\frac{2m(1/(R_0+h)-g/V^2)\cos \mathrm{\&gamma;}}{\mathrm{\&rho;S}}---\left(51\right)$

从而可得攻角和倾侧角分别满足， 

$\begin{array}{c}{C}_L(\mathrm{\&alpha;},\mathrm{Ma})={[{(C_{N0}+C_{N1})}^2+C_Y^2]}^{1/2}\\ \mathrm{\&sigma;}=\arcsin \frac{C_Y}{{[{(C_{N0}+C_{N1})}^2+C_Y^2]}^{1/2}}\end{array}---\left(52\right)$

其中：α为攻角；σ为速度滚转角；Ma为马赫数。 

（7）最优初始下滑角度解析解 

在上面的计算中，γ₀为任意给定的初值。然而由弹道阻尼理论知，合理的初始下滑角γ₀将会使得整个滑翔段的控制规律更加的平稳，具体推导如下。 

由拟平衡滑翔条件知， 

$\frac{1}{2}\mathrm{\&rho;}{V}^{2}S{C}_{N0}+m\frac{V^2}{R_0+h}\cos \mathrm{\&gamma;}-m\frac{\mathrm{\&mu;}}{{(R_0+h)}^2}\cos \mathrm{\&gamma;}=0---\left(53\right)$

其中：C_N0为拟平衡滑翔对应的纵向升力系数。假设当γ＝γ^*时，飞行器能够保持平衡滑翔状态，对上式求全微分可得， 

$\begin{array}{c}\left(\frac{1}{2}{V}^2{\mathrm{SC}}_{N0}\frac{\&PartialD;\mathrm{\&rho;}}{\&PartialD;h}\right)\frac{\mathrm{dh}}{\mathrm{dt}}+\left(\mathrm{\&rho;}{\mathrm{VSC}}_{N0}\right)\frac{\mathrm{dV}}{\mathrm{dt}}+\left(2m\frac{V}{R_0+h}\cos \mathrm{\&gamma;}\right)\frac{\mathrm{dV}}{\mathrm{dt}}\\ +(-m\frac{V^2}{{(R_0+h)}^2}\cos \mathrm{\&gamma;})\frac{\mathrm{dh}}{\mathrm{dt}}+(-m\frac{V^2}{R_0+h}\sin \mathrm{\&gamma;})\frac{\mathrm{d\&gamma;}}{\mathrm{dt}}\\ +\left(2m\frac{\mathrm{\&mu;}}{{(R_0+h)}^3}\cos \mathrm{\&gamma;}\right)\frac{\mathrm{dh}}{\mathrm{dt}}+\left(m\frac{\mathrm{\&mu;}}{{(R_0+h)}^2}\sin \mathrm{\&gamma;}\right)\frac{\mathrm{d\&gamma;}}{\mathrm{dt}}=0\end{array}---\left(54\right)$

将式(18)带入上式可得， 

${\mathrm{\&gamma;}}^{*}=\frac{\frac{D}{m}}{\frac{L\cos \mathrm{\&sigma;}{V}^2{\mathrm{\&beta;}}_h}{2\mathrm{mg}}-\frac{L\cos \mathrm{\&sigma;}}{m}-\frac{V^4}{2\mathrm{\&mu;}}}---\left(55\right)$

其中：γ^*为拟平衡滑翔最优下滑弹道倾角；Lcosσ为升力纵向分量。由于(Lcosσ/m)V²β_h/(2g)>>(Lcosσ/m)[1+mV⁴/(2μLcosσ)]，从而上式可进一步简化为， 

${\mathrm{\&gamma;}}^{*}=\frac{g}{(V^2/2){\mathrm{\&beta;}}_r\cos \mathrm{\&sigma;}}\left(\frac{1}{K_N}\right)---\left(56\right)$

其中：K_N为纵向升阻比。式(27)、式(32)、式(35)、式(42)、式(50)、式(52)和式(56)共同构成了滑翔段弹道的解析解。 

实施例 

为了检验解析求解算法的精度，选用CAV作为计算模型，进行数值仿真效验。图4至图10给出了解析弹道精度的验证结果。由图4可以看出，纵向弹道解析解的经度非常高，不管是分段多还是少，得到的结果均与数值积分完全吻合，这是由于纵向弹道求解过程中没有设定任何假设，是精确的几何弹道。同理，图5中的弹道倾角也是和数值积分完全吻合。图6可以看出，在分5段的时候解析求解的速度与数值积分得到的数度曲线还是存在微小的 差异，这是由于在实际飞行过程中不可能完全做到等纵向升阻比。由于攻角和倾侧角的精度与高度精度及速度精度相关，因此图7和图8的攻角和倾侧角曲线与速度类似（高度精度高，忽略不计）。由图9可以看出，横向弹道的精度相对较低，在分5段的时候横向弹道有明显的误差，在13段的时候误差变得较小，20段的时候精度已经较高。由图10可以看出，弹道偏角的误差较小，这说明在纬度较小情况下，纬度对弹道偏角影响较小的假设是对的。图11给出了不同初始弹道倾角情况下的解析求解弹道的对比图，由图可以看出，γ₀越大，则射程越远，在H-V走廊中的位置越靠上，但同时也会使得控制变量的震荡越厉害。因此，为了使规划出来的弹道控制上更易于实现，一般选取γ₀＝γ^*，此时对应的控制震荡最小（但还是会有震荡，因为实际的γ^*是随高度和速度变化的）。图12给出了基于该解析解的弹道曲线族。由图12(a)和图12(b)可以看出，以升力系数坡度率调整量和升力系数横向分量为控制变量的滑翔段规划弹道对应的攻角和倾侧角均较为平滑；由图12(c)和图12(d)可以看出，该方案的规划弹道纵程和横程覆盖范围均较大；由图12(e)可以看出，速度变化也较为光滑；由图12(f)、图12(g)和图12(h)可以看出，过载、热流密度和动压均满足要求；由图12(i)可以看出，基于该解析解的弹道族在约束走廊中的位置几乎重合。 

Claims (1)

1.一种高超声速飞行器滑翔段弹道解析求解方法，包括以下几个步骤：

步骤1：定义当地地心坐标系；

当地地心坐标系S_e1的坐标原点位于地心，x_e1由地心指向飞行器；y_e1在由x_e1与V的平面内，且与V同一侧，z_e1由右手定则确定；

地心赤道旋转坐标系S_e的坐标原点也为地心，z_e由地心北极，x_e处在赤道平面内指向本初子午线，y_e由右手定则确定；

（1）获取转换矩阵；

当地地心坐标系S_e1与地心赤道旋转坐标系S_e之间的转换矩阵为：

$\begin{array}{c}{L}_{e_{1}e}=L_z\left({\mathrm{\&omega;}}_1\right)L_x\left(i_1\right)L_z\left({\mathrm{\&Omega;}}_1\right)\\ =\left[\begin{array}{ccc}\left(\begin{array}{c}\cos {\mathrm{\&omega;}}_1\cos i_1\cos {\mathrm{\&Omega;}}_1\\ -\sin {\mathrm{\&omega;}}_1\sin {\mathrm{\&Omega;}}_1\end{array}\right)& \left(\begin{array}{c}\cos {\mathrm{\&omega;}}_1\cos i_1\sin {\mathrm{\&Omega;}}_1\\ +\sin {\mathrm{\&omega;}}_1\cos {\mathrm{\&Omega;}}_1\end{array}\right)& -\cos {\mathrm{\&omega;}}_1\sin i_1\\ \left(\begin{array}{c}-\sin {\mathrm{\&omega;}}_1\cos i_1\cos {\mathrm{\&Omega;}}_1\\ -\cos {\mathrm{\&omega;}}_1\sin {\mathrm{\&Omega;}}_1\end{array}\right)& \left(\begin{array}{c}-\sin {\mathrm{\&omega;}}_1\cos i_1\sin {\mathrm{\&Omega;}}_1\\ +\cos {\mathrm{\&omega;}}_1\cos {\mathrm{\&Omega;}}_1\end{array}\right)& \sin {\mathrm{\&omega;}}_1\sin i_1\\ \sin i_1\cos {\mathrm{\&Omega;}}_1& \sin i_1\sin {\mathrm{\&Omega;}}_1& \cos i_1\end{array}\right]\end{array}----\left(1\right)$

其中：

表示当地地心坐标系与地心赤道旋转坐标系之间的转换矩阵；L_z(ω₁)表示绕z轴旋转ω₁对应的转换矩阵；L_x(i₁)表示绕x轴旋转i₁对应的转换矩阵；L_z(Ω₁)表示绕z轴旋转Ω₁对应的转换矩阵；ω₁为近地点幅角，Ω₁为升交点赤经，i₁为轨道倾角，具体为：

i₁＝arccos(cosφsinψ)               (2)

${\mathrm{\&Omega;}}_1=\left\{\begin{array}{cc}\arccos \left(\frac{\cos \mathrm{\&theta;}\cos \mathrm{\&psi;}+\sin \mathrm{\&phi;}\sin \mathrm{\&psi;}\sin \mathrm{\&theta;}}{\sqrt{{\cos }^2\mathrm{\&psi;}+{\sin }^2\mathrm{\&phi;}{\sin }^2\mathrm{\&psi;}}}\right)& \frac{\sin \mathrm{\&theta;}\cos \mathrm{\&psi;}-\sin \mathrm{\&phi;}\sin \mathrm{\&psi;}\cos \mathrm{\&theta;}}{\sqrt{{\cos }^2\mathrm{\&psi;}{+\sin }^2\mathrm{\&phi;}{\sin }^2\mathrm{\&psi;}}}\&GreaterEqual;0\\ 2\mathrm{\&pi;}-\arccos \left(\frac{\cos \mathrm{\&theta;}\cos \mathrm{\&psi;}+\sin \mathrm{\&phi;}\sin \mathrm{\&psi;}\sin \mathrm{\&theta;}}{\sqrt{{\cos }^2\mathrm{\&psi;}+{\sin }^2\mathrm{\&phi;}{\sin }^2\mathrm{\&psi;}}}\right)& \frac{\sin \mathrm{\&theta;}\cos \mathrm{\&psi;}-\sin \mathrm{\&phi;}\sin \mathrm{\&psi;}\cos \mathrm{\&theta;}}{\sqrt{{\cos }^2\mathrm{\&psi;}+{\sin }^2\mathrm{\&phi;}{\sin }^2\mathrm{\&psi;}}}<0\end{array}\right.---\left(3\right)$

${\mathrm{\&omega;}}_1=\left\{\begin{array}{cc}\arccos \left(\frac{\cos \mathrm{\&phi;}\cos \mathrm{\&psi;}}{\sqrt{{\cos }^2\mathrm{\&psi;}+{\sin }^2\mathrm{\&phi;}{\sin }^2\mathrm{\&psi;}}}\right)& \mathrm{judge}>0\\ 2\mathrm{\&pi;}-\arccos \left(\frac{\cos \mathrm{\&phi;}\cos \mathrm{\&psi;}}{\sqrt{{\cos }^2\mathrm{\&psi;}+{\sin }^2\mathrm{\&phi;}{\sin }^2\mathrm{\&psi;}}}\right)& \mathrm{judge}\&le;0\end{array}\right.---\left(4\right)$

θ为在地心赤道旋转坐标系下的当地经度；φ为地心赤道旋转坐标系下的纬度；ψ为地心赤道旋转坐标系下的行向角；judge为判断近地点幅角正负号的特征量；

（2）坐标转换

设在当地地心坐标系中存在一点[θ₁,φ₁]，并且该点的弹道偏角为ψ₁，假设该点在地心赤道旋转坐标系中的坐标为[θ,φ]，弹道偏角为ψ，则它们满足，

$\left[\begin{array}{c}\cos \mathrm{\&phi;}\cos \mathrm{\&theta;}\\ \cos \mathrm{\&phi;}\sin \mathrm{\&theta;}\\ \sin \mathrm{\&phi;}\end{array}\right]={\left(L_{e_{1}e}\right)}^T\left[\begin{array}{c}\cos {\mathrm{\&phi;}}_1\cos {\mathrm{\&theta;}}_1\\ \cos {\mathrm{\&phi;}}_1\sin {\mathrm{\&theta;}}_1\\ \sin {\mathrm{\&phi;}}_1\end{array}\right]---\left(5\right)$

从而，获得在地心赤道旋转坐标系中的位置为，

$\begin{array}{c}\mathrm{\&phi;}=\arcsin (l_{11}\cos {\mathrm{\&phi;}}_1\cos {\mathrm{\&theta;}}_1+l_{21}\cos {\mathrm{\&phi;}}_1\sin {\mathrm{\&theta;}}_1+l_{31}\sin {\mathrm{\&phi;}}_1)\\ \begin{array}{}\end{array},\mathrm{\&theta;}=\left\{\begin{array}{cc}\arcsin \left(\frac{l_{12}\cos {\mathrm{\&phi;}}_1\cos {\mathrm{\&theta;}}_1+l_{22}\cos {\mathrm{\&phi;}}_1\sin {\mathrm{\&theta;}}_1+l_{32}\sin {\mathrm{\&phi;}}_1}{\cos \mathrm{\&phi;}}\right)& k_{\mathrm{\&theta;}}>0\\ \arcsin \left(\frac{l_{12}\cos {\mathrm{\&phi;}}_1\cos {\mathrm{\&theta;}}_1+l_{22}\cos {\mathrm{\&phi;}}_1\sin {\mathrm{\&theta;}}_1+l_{32}\sin {\mathrm{\&phi;}}_1}{\cos \mathrm{\&phi;}}\right)+\mathrm{\&pi;}& k_{\mathrm{\&theta;}}\&le;0\end{array}\right.\end{array}---\left(6\right)$

其中：

l₁₁＝cosω₁cosi₁cosΩ₁-sinω₁sinΩ₁

l₂₁＝-sinω₁cosi₁cosΩ₁-cosω₁sinΩ₁

l₃₁＝sini₁cosΩ₁

l₁₂＝cosω₁cosi₁sinΩ₁+sinω₁cosΩ₁

l₂₂＝-sinω₁cosi₁sinΩ₁+cosω₁cosΩ₁

l₃₂＝sini₁sinΩ₁

k_θ＝l₁₁cosφ₁cosθ₁+l₂₁cosφ₁sinθ₁+l₃₁sinφ₁

单位速度矢量定义如下，

$\left[\begin{array}{c}{v}_x\\ v_y\\ v_z\end{array}\right]={\left(L_{e_{1}e}\right)}^T\left[[\begin{array}{c}-\sin {\mathrm{\&psi;}}_1\sin {\mathrm{\&theta;}}_1-\cos {\mathrm{\&psi;}}_1\sin {\mathrm{\&phi;}}_1\cos {\mathrm{\&theta;}}_1\\ \sin {\mathrm{\&psi;}}_1\cos {\mathrm{\&theta;}}_1-\cos {\mathrm{\&psi;}}_1\sin {\mathrm{\&phi;}}_1\sin {\mathrm{\&theta;}}_1\\ \cos {\mathrm{\&psi;}}_1\cos {\mathrm{\&phi;}}_1\end{array}\right]---\left(7\right)$

另外，在地心赤道旋转坐标系中，

$\left[\begin{array}{c}-\sin \mathrm{\&psi;}\sin \mathrm{\&theta;}-\cos \mathrm{\&psi;}\sin \mathrm{\&phi;}\cos \mathrm{\&theta;}\\ \sin \mathrm{\&psi;}\cos \mathrm{\&theta;}-\cos \mathrm{\&psi;}\sin \mathrm{\&phi;}\sin \mathrm{\&theta;}\\ \cos \mathrm{\&psi;}\cos \mathrm{\&phi;}\end{array}\right]=\left[\begin{array}{c}{v}_x\\ v_y\\ v_z\end{array}\right]$

由式(7)得，

$\begin{array}{c}\cos \mathrm{\&psi;}=\frac{v_z}{\cos \mathrm{\&phi;}}\\ \sin \mathrm{\&psi;}=\left\{\begin{array}{cc}\frac{v_x+\cos \mathrm{\&psi;}\sin \mathrm{\&phi;}\cos \mathrm{\&theta;}}{-\sin \mathrm{\&theta;}}& 0\&Element;\{0,\mathrm{\&pi;}\}\\ \frac{v_y+\cos \mathrm{\&psi;}\sin \mathrm{\&phi;}\sin \mathrm{\&theta;}}{\cos \mathrm{\&theta;}}& 0\&Element;(0,\mathrm{\&pi;})\&cup;(\mathrm{\&pi;},2\mathrm{\&pi;})\end{array}\right.\end{array}---\left(8\right)$

从而可得，

$\mathrm{\&psi;}=\left\{\begin{array}{cc}\arccos \left(\frac{v_z}{\cos \mathrm{\&phi;}}\right)& \sin \mathrm{\&psi;}\&GreaterEqual;0\\ 2\mathrm{\&pi;}-\arccos \left(\frac{v_z}{\cos \mathrm{\&phi;}}\right)& \sin \mathrm{\&psi;}<0\end{array}\right.---\left(9\right)$

其中：v_y为单位速度矢量y轴分量；v_z为单位速度矢量z轴分量；

步骤2：获取在当地地心坐标系下的运动方程

无动力滑翔飞行器单位质量机械能表述为e＝μ/r-V²/2，忽略地球自转，在当地地心坐标系中三自由度运动模型可表述为，

$\begin{array}{c}\stackrel{\&CenterDot;}{h}=V\sin \mathrm{\&gamma;}\\ {\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_1=\frac{V\cos \mathrm{\&gamma;}\sin {\mathrm{\&psi;}}_1}{(h+R_0)\cos {\mathrm{\&phi;}}_1}\\ {\stackrel{\&CenterDot;}{\mathrm{\&phi;}}}_1=\frac{V\cos \mathrm{\&gamma;}\cos {\mathrm{\&psi;}}_1}{h+R_0}\\ \stackrel{\&CenterDot;}{e}=\frac{\mathrm{DV}}{m}\\ \stackrel{\&CenterDot;}{\mathrm{\&gamma;}}=\frac{1}{V}[\frac{L\cos \mathrm{\&sigma;}}{m}+(\frac{V^2}{h+R_0}-g)\cos \mathrm{\&gamma;}]\\ {\stackrel{\&CenterDot;}{\mathrm{\&psi;}}}_1=\frac{1}{V}[\frac{L\sin \mathrm{\&sigma;}}{m\cos \mathrm{\&gamma;}}+\frac{V^2}{h+R_0}\cos \mathrm{\&gamma;}\sin {\mathrm{\&psi;}}_1\tan {\mathrm{\&phi;}}_1]\\ \stackrel{\&CenterDot;}{s}=\frac{V\cos \mathrm{\&gamma;}}{h+R_0}\end{array}---\left(10\right)$

其中：R₀＝6378km，θ₁和φ₁为当地地心坐标系下的经度和纬度；V为飞行器与地球的相对速度；ψ₁为当地地心坐标系下的弹道偏角；γ为地心赤道旋转坐标系下的弹道倾角；h为地心赤道旋转坐标系下的高度；s为滑翔段的飞行路程；

和分别为高度、当地地心坐标系下的经度、当地地心坐标系下的纬度、负比能量、弹道倾角和当地地心坐标系下的弹道偏角及射程对时间的一阶导数；D和L分别为以过载形式表示的升力和阻力，即D＝0.5ρV²SC_D，L＝0.5ρV²SC_L，其中ρ为大气密度，S为气动参考面积，m为飞行器的质量；C_L和C_D分别为升力系数和阻力系数；g为当地重力加速度；σ为速度滚转角；

由于φ₁≈0、γ≈0、ψ₁≈π/2，并假设0.5ρV²SC_N0＝-m(V²/r-g)cosγ且C_Lcosσ＝C_N0+C_N1，其中C_N0为纵向平衡滑翔升力系数、C_N1为升力系数坡度率调整量，式(10)可进一步化为：

$\begin{array}{c}\begin{array}{cc}\stackrel{\&CenterDot;}{h}=V\sin & \stackrel{\&CenterDot;}{e}=\frac{\mathrm{DV}}{m}\\ {\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_1=\frac{V}{h+R_0}& \stackrel{\&CenterDot;}{\mathrm{\&gamma;}}=\frac{\mathrm{\&rho;}{\mathrm{VSC}}_{N1}}{2m}\\ {\stackrel{\&CenterDot;}{\mathrm{\&phi;}}}_1=\frac{V(\mathrm{\&pi;}/2-{\mathrm{\&psi;}}_1)}{h+R_0}& {\stackrel{\&CenterDot;}{\mathrm{\&psi;}}}_1=\frac{\mathrm{\&rho;}{\mathrm{VSC}}_Y}{2m}+\frac{V}{h+R_0}{\mathrm{\&phi;}}_1\end{array}\\ \stackrel{\&CenterDot;}{s}=\frac{V\cos \mathrm{\&gamma;}}{h+R_0}\end{array}---\left(11\right)$

其中：C_N1为升力系数坡度率调整量，C_Y为升力系数横向调整量，满足C_Y＝C_Lsinσ；

取ψ₂＝π/2-ψ₁，并对上式进行如下的变换，

$\frac{\mathrm{dh}}{\sin \mathrm{\&gamma;}}=\mathrm{Vdt}---\left(12\right)$

(h+R₀)dθ₁＝Vdt              (13)

$\frac{(h+R_0)d{\mathrm{\&phi;}}_1}{{\mathrm{\&psi;}}_2}=\mathrm{Vdt}---\left(14\right)$

$\frac{2\mathrm{md\&gamma;}}{\mathrm{\&rho;}{\mathrm{SC}}_{N1}}=\mathrm{Vdt}---\left(15\right)$

$\frac{d{\mathrm{\&psi;}}_2}{-\frac{\mathrm{\&rho;}{\mathrm{SC}}_Y}{2m}-\frac{1}{h+R_0}{\mathrm{\&phi;}}_1}=\mathrm{Vdt}---\left(16\right)$

$\frac{\mathrm{mde}}{D}=\mathrm{Vdt}---\left(17\right)$

(h+R₀)ds＝Vdt              (18)

其中：dh、dθ₁、dφ₁、de、dγ、dψ₁和ds分别为高度、当地地心坐标系下的经度、当地地心坐标系下的纬度、负比能量、弹道倾角和当地地心坐标系下的弹道偏角及射程的微分，dt为时间的微分；

步骤3：获取滑翔弹道的解析解

设起滑点坐标为(θ₀,φ₀)，起滑高度为h₀，起滑弹道倾角为γ₀，终端高度为h_f，纵向升力系数余量为C_N1，横向升力系数为C_Y；利用式(1)可获得该点处当地地心坐标系与地心赤道旋转坐标系之间的转换关系，并且该点在当地地心坐标系中满足θ₁₀＝0、φ₁₀＝0以及ψ₁₀＝π/2；则在当地地心坐标系下滑翔段弹道的解析解如下；

（1）滑翔弹道倾角解析解

由式(12)和式(15)可得，

$\frac{\mathrm{dh}}{\sin \mathrm{\&gamma;}}=\frac{2\mathrm{md\&gamma;}}{\mathrm{\&rho;}{\mathrm{SC}}_{N1}}---\left(19\right)$

对式(19)进行积分可得，

${\mathrm{\&gamma;}}_f-\cos {\mathrm{\&gamma;}}_0=\frac{{\mathrm{SC}}_{N1}}{2m{\mathrm{\&beta;}}_h}({\mathrm{\&rho;}}_f-{\mathrm{\&rho;}}_0)---\left(20\right)$

其中：ρ_f为预测终点大气密度；γ_f为预测终点弹道倾角；β_h为指数大气模型常数；对式(20)进一步处理可得，

$\frac{{\mathrm{SC}}_{N1}}{2m{\mathrm{\&beta;}}_h}{\mathrm{\&rho;}}_f-\cos {\mathrm{\&gamma;}}_f=\frac{{\mathrm{SC}}_{N1}}{2m{\mathrm{\&beta;}}_h}{\mathrm{\&rho;}}_0-\cos {\mathrm{\&gamma;}}_0---\left(21\right)$

从而得到常数，

$K^{*}=\frac{{\mathrm{SC}}_{N1}{\mathrm{\&rho;}}_0}{2m{\mathrm{\&beta;}}_h}-\cos {\mathrm{\&gamma;}}_0---\left(22\right)$

上述的K^*即为决定滑翔段纵向弹道弯曲程度的常数；

（2）飞行路程解析解

由式(15)和式(18)可得：

$(h+R_0)\mathrm{ds}=\frac{2\mathrm{md\&gamma;}}{\mathrm{\&rho;}{\mathrm{SC}}_{N1}}---\left(23\right)$

将式(21)和式(22)带入式(23)可得：

${\mathrm{\&beta;}}_h(h+R_0)\mathrm{ds}=\frac{1}{\cos \mathrm{\&gamma;}+K^{*}}\mathrm{d\&gamma;}---\left(24\right)$

设

为地心平均距离；对式(24)进行积分可得：

${\mathrm{\&beta;}}_h\stackrel{\&OverBar;}{R}{s}_f=\left\{\begin{array}{cc}\frac{1}{1+K^{*}}\sqrt{\frac{1+K^{*}}{1-K^{*}}}\left(\begin{array}{c}\ln \left|\frac{\tan \frac{{\mathrm{\&gamma;}}_f}{2}+\sqrt{\frac{1+K^{*}}{1-K^{*}}}}{\tan \frac{{\mathrm{\&gamma;}}_f}{2}-\sqrt{\frac{1+K^{*}}{1-K^{*}}}}\right|\\ -\ln |\frac{\tan \frac{{\mathrm{\&gamma;}}_0}{2}+\sqrt{\frac{1+K^{*}}{1-K^{*}}}}{\tan \frac{{\mathrm{\&gamma;}}_0}{2}-\sqrt{\frac{1+K^{*}}{1-K^{*}}}}\end{array}\right)& \left|K^{*}\right|<1\\ \frac{2}{(1+K^{*})}\sqrt{\frac{1+K^{*}}{K^{*}-1}}\left(\begin{array}{c}a\tan \left(\sqrt{\frac{K^{*}-1}{K^{*}+1}}\tan \frac{{\mathrm{\&gamma;}}_f}{2}\right)\\ -a\tan \left(\sqrt{\frac{K^{*}-1}{K^{*}+1}}\tan \frac{{\mathrm{\&gamma;}}_0}{2}\right)\end{array}\right)& \left|K^{*}\right|>1\\ \frac{2}{K^{*}-1}(\cot \frac{{\mathrm{\&gamma;}}_f}{2}-\cot \frac{{\mathrm{\&gamma;}}_0}{2})& \left|K^{*}\right|=1\end{array}\right.---\left(25\right)$

其中：s_f为预测终点航程；

（3）滑翔经度解析解

由式(13)和式(15)可得，

$(h+R_0)d{\mathrm{\&theta;}}_1=\frac{2\mathrm{md\&gamma;}}{\mathrm{\&rho;}{\mathrm{SC}}_{N1}}---\left(26\right)$

式(26)与式(23)比较可得，

dθ₁＝ds             (27)

从而可得，

θ_1f＝s_f             (28)

其中：θ_1f为当地地心坐标系下的预测终点经度；

（4）滑翔纬度与弹道偏角解析解

由式(14)、式(15)和式(16)可得，

$\begin{array}{c}\frac{d{\mathrm{\&phi;}}_1}{\mathrm{d\&gamma;}}=\frac{2m}{(h+R_0)\mathrm{\&rho;}{\mathrm{SC}}_{N1}}{\mathrm{\&psi;}}_2\\ \frac{d{\mathrm{\&psi;}}_2}{\mathrm{d\&gamma;}}=-\frac{C_Y}{C_{N1}}{-}_{}\frac{2m}{(h+R_0){\mathrm{\&rho;SC}}_{N1}}{\mathrm{\&phi;}}_1\end{array}---\left(29\right)$

将式(21)和式(22)带入式(29)，可得，

$\begin{array}{c}\frac{d{\mathrm{\&phi;}}_1}{\mathrm{d\&gamma;}}=\frac{{\mathrm{\&psi;}}_2}{\stackrel{\&OverBar;}{R}{\mathrm{\&beta;}}_h\left(K^{*}+\cos \mathrm{\&gamma;}\right)}\\ \frac{\mathrm{d\&psi;}{2}_{}}{\mathrm{d\&gamma;}}=-\frac{C_Y}{C_{N1}}-\frac{{\mathrm{\&phi;}}_1}{\stackrel{\&OverBar;}{R}{\mathrm{\&beta;}}_h(K^{*}+\cos \mathrm{\&gamma;})}\end{array}---\left(30\right)$

取，

$\begin{array}{c}{K}_1=\frac{1}{{\mathrm{\&gamma;}}_f-{\mathrm{\&gamma;}}_0}{\&Integral;}_{{\mathrm{\&gamma;}}_0}^{{\mathrm{\&gamma;}}_f}\frac{1}{\stackrel{\&OverBar;}{R}{\mathrm{\&beta;}}_h(K^{*}+\cos \mathrm{\&gamma;})}\mathrm{d\&gamma;}\\ =\frac{{\mathrm{\&theta;}}_{1f}}{{\mathrm{\&gamma;}}_f-{\mathrm{\&gamma;}}_0}\end{array}---\left(31\right)$

则，式（30）可化为，

$\begin{array}{c}\frac{d{\mathrm{\&phi;}}_1}{\mathrm{d\&gamma;}}=K_1{\mathrm{\&psi;}}_2\\ \frac{d{\mathrm{\&psi;}}_2}{\mathrm{d\&gamma;}}=-\frac{C_Y}{C_{N1}}-K_1{\mathrm{\&phi;}}_1\end{array}---\left(32\right)$

由于φ₁₀＝0，ψ₂₀＝0，对上式进行拉式变换可得，

$\begin{array}{c}s{\mathrm{\&phi;}}_1\left(s\right)=K_1{\mathrm{\&psi;}}_2\left(s\right)\\ s{\mathrm{\&psi;}}_2\left(s\right)=-\frac{C_Y}{C_{N1}}\frac{1}{s}-K_1{\mathrm{\&phi;}}_1\left(s\right)\end{array}---\left(33\right)$

由上式可得，

$\begin{array}{c}{\mathrm{\&psi;}}_2\left(s\right)=\frac{-\frac{C_Y}{C_{N1}}}{(s^2+K_1^2)}\\ {\mathrm{\&phi;}}_1\left(s\right)=\frac{-\frac{C_Y}{C_{N1}}{K}_1}{s(s^2+K_1^2)}\end{array}---\left(34\right)$

对式(34)进行拉式反变换可得，

$\begin{array}{c}{\mathrm{\&psi;}}_{2f}=-\frac{C_Y}{C_{N1}}\frac{1}{K_1}\sin {\mathrm{\&theta;}}_{1f}\\ {\mathrm{\&phi;}}_{1f}=-\frac{C_Y}{C_{N1}}\frac{1}{K_1}(1-\cos {\mathrm{\&theta;}}_{1f})\end{array}---\left(35\right)$

其中：φ_1f为当地地心坐标系下的预测终点纬度；由上式还可得当地地心坐标系下的预测终点弹道偏角ψ_1f为ψ_1f＝π/2-ψ_2f；利用式(35)和式(1)给出的坐标转换关系，就可以获得在地心赤道旋转坐标系下预测终点经度θ_f预测终点纬度φ_f和预测终点弹道偏角ψ_f；

（5）滑翔速度解析解

假定在滑翔过程中，保持定常纵向升阻比K_N，K_N＝C_Lcosσ/C_D，则滑翔阻力为，

$D=\frac{L\cos \mathrm{\&sigma;}}{K_N}---\left(36\right)$

由于Lcosσ＝0.5ρV²S(C_N0+C_N1)，忽略C_N1后可得，

$D=\frac{m[g-V^2/(R_0+h)]\cos \mathrm{\&gamma;}}{K_N}+\frac{0.5\mathrm{\&rho;}{V}^2{\mathrm{SC}}_{N1}}{K_N}---\left(37\right)$

其中：D为阻力；m为飞行器质量；V为飞行器相对地球的速度；C_N1为升力系数坡度率调整量；g为当地重力加速度，写成如下形式，

$g=\frac{\mathrm{\&mu;}}{{(R_0+h)}^2}$

其中：μ为由地球引力常数；

式(37)中，阻力可分为两个部分，即保持平衡滑翔的升力产生的阻力以及用于纵向转弯的升力产生的阻力，忽略第二部分的影响，第一部分阻力满足，

$\frac{\mathrm{de}}{d{\mathrm{\&theta;}}_1}=\frac{-\frac{\mathrm{\&mu;}}{R_0+h}+2e}{K_N}---\left(38\right)$

由于h相对R₀为小量，采用它的算术平均值代替，则上式积分可得，

$\ln \left(\frac{2e-\mathrm{\&mu;}/\stackrel{\&OverBar;}{R}}{2e_0-\mathrm{\&mu;}/\stackrel{\&OverBar;}{R}}\right)=\frac{2{\mathrm{\&theta;}}_{1f}}{K_N}---\left(39\right)$

其中：e₀为初始负比能量。从而可得，

$2e_f=\mathrm{\&mu;}/\stackrel{\&OverBar;}{R}+(2e_0-\mathrm{\&mu;}/\stackrel{\&OverBar;}{R})\exp \left(\frac{2{\mathrm{\&theta;}}_{1f}}{K}\right)---\left(40\right)$

其中：e_f为预测终点负比能量；对于另一部分用于转弯的升力，其消耗的机械能为，

$\mathrm{\&Delta;e}={\&Integral;}_{s_0}^s\frac{V\stackrel{\&CenterDot;}{\mathrm{\&gamma;}}}{K_N}\mathrm{rds}={\&Integral;}_t^{t_f}\frac{V^2\stackrel{\&CenterDot;}{\mathrm{\&gamma;}}}{K_N}\mathrm{dt}---\left(41\right)$

其中：为弹道倾角的导数；

取 ${\stackrel{\&OverBar;}{V}}^2=(V_0^2+V_f^2+V_0{V}_f)/3,$ 则可得，

$\mathrm{\&Delta;e}=\frac{{\stackrel{\&OverBar;}{V}}^2}{K_N}({\mathrm{\&gamma;}}_f-{\mathrm{\&gamma;}}_0)---\left(42\right)$

其中：Δe为负比能量修正量，综上可得，终端速度预测为：

$V_f=\sqrt{\mathrm{\&mu;}/\stackrel{\&OverBar;}{R}-(2e_0-\mathrm{\&mu;}/\stackrel{\&OverBar;}{R})\exp (2{\mathrm{\&theta;}}_{1f}/K_N)-\frac{2{\stackrel{\&OverBar;}{V}}^2}{K_N}({\mathrm{\&gamma;}}_f-{\mathrm{\&gamma;}}_0)}---\left(43\right)$

其中：V_f为预测终点速度；

为距离地心的平均距离；K_N为升阻比；

在式(43)的计算中利用了平均速度，需要通过迭代来不断修正

同时，利用迭代修正对K_N的估计；

（6）滑翔攻角和倾侧角解析解

式(21)、式(28)、式(35)、式(43)分别给出了预测终点弹道倾角γ_f、经度θ_f、纬度φ_f、弹道偏角ψ_f、速度V_f与高度h_f之间的关系；利用平衡滑翔关系可知，

$C_{N0}=-\frac{2m(1/(R_0+h)-g/V^2)\cos \mathrm{\&gamma;}}{\mathrm{\&rho;S}}---\left(44\right)$

从而可得攻角和倾侧角分别满足，

$\begin{array}{c}{C}_L(\mathrm{\&alpha;},\mathrm{Ma})={[{(C_{N0}+C_{N1})}^2+C_Y^2]}^{1/2}\\ \mathrm{\&sigma;}=\arcsin \frac{C_Y}{{[{(C_{N0}+C_{N1})}^2+C_Y^2]}^{1/2}}\end{array}---\left(45\right)$

其中：α为攻角；σ为速度滚转角；Ma为马赫数；

（7）最优初始下滑角度解析解

具体推导如下；

由拟平衡滑翔条件知，

$\frac{1}{2}\mathrm{\&rho;}{V}^{2}S{C}_{N0}+m\frac{V^2}{R_0+h}\cos \mathrm{\&gamma;}-m\frac{\mathrm{\&mu;}}{{(R_0+h)}^2}\cos \mathrm{\&gamma;}=0---\left(46\right)$

其中：C_N0为拟平衡滑翔对应的纵向升力系数；假设当γ＝γ^*时，飞行器能够保持平衡滑翔状态，对上式求全微分可得，

$\begin{array}{c}\left(\frac{1}{2}{V}^2{\mathrm{SC}}_{N0}\frac{\&PartialD;\mathrm{\&rho;}}{\&PartialD;h}\right)\frac{\mathrm{dh}}{\mathrm{dt}}+\left(\mathrm{\&rho;}{\mathrm{VSC}}_{N0}\right)\frac{\mathrm{dV}}{\mathrm{dt}}+\left(2m\frac{V}{R_0+h}\cos \mathrm{\&gamma;}\right)\frac{\mathrm{dV}}{\mathrm{dt}}\\ +(-m\frac{V^2}{{(R_0+h)}^2}\cos \mathrm{\&gamma;})\frac{\mathrm{dh}}{\mathrm{dt}}+(-m\frac{V^2}{R_0+h}\sin \mathrm{\&gamma;})\frac{\mathrm{d\&gamma;}}{\mathrm{dt}}\\ +\left(2m\frac{\mathrm{\&mu;}}{{(R_0+h)}^3}\cos \mathrm{\&gamma;}\right)\frac{\mathrm{dh}}{\mathrm{dt}}+\left(m\frac{\mathrm{\&mu;}}{{(R_0+h)}^2}\sin \mathrm{\&gamma;}\right)\frac{\mathrm{d\&gamma;}}{\mathrm{dt}}=0\end{array}---\left(47\right)$

将式(11)带入上式可得，

${\mathrm{\&gamma;}}^{*}=\frac{\frac{D}{m}}{\frac{L\cos \mathrm{\&sigma;}{V}^2{\mathrm{\&beta;}}_h}{2\mathrm{mg}}-\frac{L\cos \mathrm{\&sigma;}}{m}-\frac{V^4}{2\mathrm{\&mu;}}}---\left(48\right)$

其中：γ^*为拟平衡滑翔最优下滑弹道倾角；Lcosσ为升力纵向分量；由于(Lcosσ/m)V²β_h(2g)>>(Lcosσm)[1+mV⁴/(2μLcosσ)]，从而上式可进一步简化为：

${\mathrm{\&gamma;}}^{*}=\frac{g}{(V^2/2){\mathrm{\&beta;}}_r\cos \mathrm{\&sigma;}}\left(\frac{1}{K_N}\right)---\left(49\right)$

其中：K_N为纵向升阻比；式(20)、式(25)、式(28)、式(35)、式(43)、式(45)和式(49)共同构成了滑翔段弹道的解析解。

Images