基于摄动理论的滑翔弹道误差传播分析方法

摘要

本发明涉及一种基于摄动理论的滑翔弹道误差传播分析方法。该方法首先以滑翔弹道再入纵平面为换极赤道平面建立换极坐标系，其次建立换极坐标系中的飞行器动力学模型及运动摄动方程，然后通过求解状态转移矩阵获得摄动方程的半解析解，最后通过坐标变化建立一般坐标系中的误差传播模型。本发明提供的方法力争弥补高超声速滑翔飞行器滑翔弹道摄动因素影响机理分析方法上的空白，为深化误差传播机理认识、建立弹道快速修正方法提供基础和方法支撑。

说明书

技术领域

本发明属于飞行器动力学建模领域，特别涉及高超声速滑翔飞行器滑翔弹道摄动因素 误差传播机理分析方法。

背景技术

飞行器弹道规划及制导控制设计通常由于物理认知不完善或为提高计算效率而基于 近似模型来开展。由近似而引入的误差主要涉及飞行器本体建模误差及地球物理环境因素 误差，统称为影响飞行状态的摄动因素。摄动因素将导致飞行器实际运动状态偏离设计状 态，是弹道规划及制导控制设计中方法误差的主要来源。摄动因素主要包括飞行器本体建 模误差和地球物理环境因素误差两类。其中，飞行器本体建模误差主要包括弹体结构误差、 发动机特性误差和气动系数误差等，地球物理环境因素误差主要包括大气模型不确定性及 扰动引力等。

削弱摄动因素影响的思路有两类，其一是通过精细化建模来削弱摄动因素本身，其二 是首先基于近似模型进行初步设计，之后进行摄动因素影响机理分析，最终建立等效补偿 方法。前者通常由于知识体系不完善、认知手段局限以及采用精细化模型将导致不可容忍 的计算量等原因而不具有理论可行性及工程应用价值。作为解决此类问题的合理手段，后 者的难点在于获得对摄动因素误差传播机理的认识，因此亟待建立相应的分析方法，以作 为寻求等效补偿机制及补偿量的基础。从现有文献来看，尚无针对高超声速滑翔飞行器建 立的摄动因素影响分析方法，使得相近领域的理论和方法难以直接应用。

发明内容

针对上述问题，本发明提出一种基于摄动理论的滑翔弹道误差传播分析方法。该方法 首先以滑翔弹道再入纵平面为换极赤道平面建立换极坐标系，其次建立换极坐标系中的飞 行器动力学模型及运动摄动方程，然后通过求解状态转移矩阵获得摄动方程的半解析解， 最后通过坐标变化建立一般坐标系中的误差传播模型。本发明提供的方法力争弥补高超声 速滑翔飞行器滑翔弹道摄动因素影响机理分析方法上的空白，为深化误差传播机理认识、 建立弹道快速修正方法提供基础和方法支撑。

本发明针对高超声速滑翔飞行器滑翔弹道的摄动因素误差传播机理问题，首次提出一 种基于摄动理论的滑翔弹道误差传播分析方法。

该方法是针对高超声速滑翔飞行器滑翔弹道建立的首个误差机理分析方法，综合利用 坐标系转换、飞行器动力学建模、摄动理论和高维常微分方程求解技术等实现弹道误差传 播模型的建立。建立的弹道误差传播模型立足于飞行器动力学方程及弹道特性，本质上是 以飞行状态偏差量为状态变量建立的动力系统，其输出用于表征飞行器对外部摄动因素的 响应特性，可广泛适用于所有以这类动力学模型描述的飞行器而不受飞行器特征参数改变 的影响。模型可独立描述初态误差和飞行过程中摄动因素的影响，可适用于单一摄动因素 和复杂摄动因素的影响分析。

本发明的误差传播模型建立步骤如下：第一步，换极坐标系的建立；第二步，换极坐 标系中飞行器动力学模型的建立；第三步，一般坐标系与换极坐标系中飞行状态量的坐标 转换；第四步，摄动方程的建立；第五步，摄动方程状态转移矩阵Φ(t_k,t)的求解；第六步，

一般坐标系中误差传播模型的建立。

本发明中，换极是极点变换的意思，换极坐标系是指基于极点变换思想以重新定义的 换极赤道面为基准建立的坐标系，而一般坐标系是指以地球赤道面为基准建立的坐标系。

具体地：本发明的技术方案包括以下步骤：

第一步，换极坐标系建立

按如下方式建立换极坐标系：

①定义一个再入大圆弧平面作为换极赤道平面：1)对目标点确定的情况，将滑翔起 点和目标点地心矢径构成的再入大圆弧平面作为换极赤道平面；2)对于目标点未确定的 情况，根据滑翔起点位置及方位角确定的再入大圆弧平面作为换极赤道面。

②基于换极赤道平面定义换极坐标系O_E为地心，轴沿滑翔起点地心矢 径方向，轴在换极赤道面内垂直于轴指向目标点方向，轴与轴、轴构成右手系。

第二步，换极坐标系中飞行器动力学模型建立

以高超声速滑翔飞行器滑翔弹道为研究对象，确定其飞行状态量为经度地心纬度航迹偏航角速度速度倾角和地心距建立以时间为自变量的动力学模型，

$\left(\begin{array}{c}\frac{d\hat{r}}{dt}=\hat{V}\sin \hat{\theta }\\ \frac{d\hat{\lambda }}{dt}=\frac{\hat{V}\cos \hat{\theta }\sin \hat{\sigma }}{\hat{r}\cos \hat{\phi }}\\ \frac{d\hat{\phi }}{dt}=\frac{\hat{V}\cos \hat{\theta }\cos \hat{\sigma }}{\hat{r}}\\ \frac{d\hat{\sigma }}{dt}=\frac{L\sin \upsilon }{\hat{V}\cos \hat{\theta }}+\frac{\hat{V}}{\hat{r}}\cos \hat{\theta }\tan \hat{\phi }\sin \hat{\sigma }-{\hat{g}}_{{\omega }_e}\frac{\cos \hat{\phi }\sin \hat{\sigma }}{\hat{V}\cos \hat{\theta }}+C_{\sigma }+{\tilde{C}}_{\sigma }\\ \frac{d\hat{\theta }}{dt}=\frac{L}{\hat{V}}\cos \upsilon +\frac{\hat{V}}{\hat{r}}\cos \hat{\theta }+\frac{{\hat{g}}_r\cos \hat{\theta }}{\hat{V}}+\frac{{\hat{g}}_{{\omega }_e}}{\hat{V}}(\cos \hat{\theta }\sin \hat{\phi }-\cos \hat{\sigma }\sin \hat{\theta }\cos \hat{\phi })+C_{\theta }+{\tilde{C}}_{\theta }\\ \frac{d\hat{V}}{dt}=-D+{\hat{g}}_r\sin \hat{\theta }+{\hat{g}}_{{\omega }_e}(\cos \hat{\sigma }\cos \hat{\theta }\cos \hat{\phi }+\sin \hat{\theta }\sin \hat{\phi })+{\tilde{C}}_V\end{array}\right)---\left(25\right)$

其中，C_σ、C_θ为哥氏加速度项，和为牵连加速度项，其表达式如下，

$\left(\begin{array}{c}{C}_{\sigma }=(2{\omega }_{ex}-2\tan \hat{\theta }({\omega }_{ey}\sin \hat{\sigma }-{\omega }_{ez}\sin \hat{\sigma }\left)\right)\\ {\tilde{C}}_{\sigma }=-\frac{\hat{r}}{\hat{V}\cos \hat{\theta }}({\omega }_{ex}{\omega }_{ey}\cos \hat{\sigma }-{\omega }_{ex}{\omega }_{ez}\sin \hat{\sigma })\\ C_{\theta }=2({\omega }_{ez}\sin \hat{\sigma }-{\omega }_{ey}\cos \hat{\sigma })\\ {\tilde{C}}_{\theta }=\frac{\hat{r}}{\hat{V}}[{\omega }_{ex}{\omega }_{ey}\sin \hat{\theta }\sin \hat{\sigma }+{\omega }_{ex}{\omega }_{ez}\sin \hat{\theta }\cos \hat{\sigma }+({\omega }_{ey}^2+{\omega }_{ez}^2)\cos \hat{\theta }]\\ {\tilde{C}}_V=\hat{r}[-{\omega }_{ex}{\omega }_{ey}soc\hat{\theta }\sin \hat{\sigma }-{\omega }_{ex}{\omega }_{ez}soc\hat{\theta }\cos \hat{\sigma }+({\omega }_{ey}^2+{\omega }_{ez}^2)\sin \hat{\theta }]\end{array}\right)---\left(26\right)$

其中，

$\left(\begin{array}{c}{\omega }_{ex}={\omega }_e(\cos \hat{\lambda }\cos \hat{\phi }{\mathrm{cos\phi }}_p{\mathrm{cosA}}_p+\sin \hat{\lambda }\cos \hat{\phi }{\mathrm{cos\phi }}_p{\mathrm{sinA}}_p+\sin \hat{\phi }{\mathrm{sin\phi }}_p)\\ {\omega }_{ey}={\omega }_e(-\sin \hat{\lambda }{\mathrm{cos\phi }}_p{\mathrm{cosA}}_p+\cos \hat{\lambda }{\mathrm{cos\phi }}_p{\mathrm{sinA}}_p)\\ {\omega }_{ez}={\omega }_e(-\cos \hat{\lambda }\sin \hat{\phi }{\mathrm{cos\phi }}_p{\mathrm{cosA}}_p-\sin \hat{\lambda }\sin \hat{\phi }{\mathrm{cos\phi }}_p{\mathrm{sinA}}_p+\cos \hat{\phi }{\mathrm{sin\phi }}_p)\end{array}\right)---\left(27\right)$

其中，ω_e为地球旋转加速度矢量，λ_p为换极后极点P的经度，φ_p为P的地心纬度，A_P为P的方位角。

第三步，一般坐标系与换极坐标系中飞行状态量坐标转换

根据换极坐标系定义，一般坐标系与换极坐标系中地心距、当地速度倾角及速度的定 义一致，

$r=\hat{r},\theta =\hat{\theta },V=\hat{V}---\left(28\right)$

定义

$\begin{array}{c}\left(\begin{array}{ccc}{\mathrm{cos\phi }}_0{\mathrm{cos\lambda }}_0& {\mathrm{cos\phi }}_0{\mathrm{sin\lambda }}_0& {\mathrm{sin\phi }}_0\\ -{\mathrm{sin\psi }}_0{\mathrm{sin\lambda }}_0-{\mathrm{cos\psi }}_0{\mathrm{sin\phi }}_0{\mathrm{cos\lambda }}_0& {\mathrm{sin\psi }}_0{\mathrm{cos\lambda }}_0-{\mathrm{cos\psi }}_0{\mathrm{sin\phi }}_0{\mathrm{sin\lambda }}_0& {\mathrm{cos\psi }}_0{\mathrm{cos\phi }}_0\\ {\mathrm{cos\psi }}_0{\mathrm{sin\lambda }}_0-{\mathrm{sin\psi }}_0{\mathrm{sin\phi }}_0{\mathrm{cos\lambda }}_0& -{\mathrm{cos\psi }}_0{\mathrm{cos\lambda }}_0-{\mathrm{sin\psi }}_0{\mathrm{sin\phi }}_0{\mathrm{sin\lambda }}_0& {\mathrm{sin\psi }}_0{\mathrm{cos\phi }}_0\end{array}\right)\\ \stackrel{\Delta }{=}\left(\begin{array}{ccc}{G}_{11}& G_{12}& G_{13}\\ G_{21}& G_{22}& G_{23}\\ G_{31}& G_{32}& G_{33}\end{array}\right).\end{array}---\left(29\right)$

$\{\begin{array}{c}{G}_{11}\cos \phi \cos \lambda +G_{12}\cos \phi \sin \lambda +G_{13}\sin \phi \stackrel{\Delta }{=}{k}_1\\ G_{21}\cos \phi \cos \lambda +G_{22}\cos \phi \sin \lambda +G_{23}\sin \phi \stackrel{\Delta }{=}{k}_2\\ G_{31}\cos \phi \cos \lambda +G_{32}\cos \phi \sin \lambda +G_{33}\sin \phi \stackrel{\Delta }{=}{k}_3\end{array}---\left(30\right)$

$\{\begin{array}{c}{G}_{11}\cos \hat{\phi }\cos \hat{\lambda }+G_{21}\cos \hat{\phi }\sin \hat{\lambda }+G_{31}\sin \hat{\phi }\stackrel{\Delta }{=}{\tilde{k}}_1\\ G_{12}\cos \hat{\phi }\cos \hat{\lambda }+G_{22}\cos \hat{\phi }\sin \hat{\lambda }+G_{32}\sin \hat{\phi }\stackrel{\Delta }{=}{\tilde{k}}_2\\ G_{13}\cos \hat{\phi }\cos \hat{\lambda }+G_{23}\cos \hat{\phi }\sin \hat{\lambda }+G_{33}\sin \hat{\phi }\stackrel{\Delta }{=}{\tilde{k}}_3\end{array}---\left(31\right)$

由λ和φ确定和的表达式为，

$\begin{array}{cc}\left(\begin{array}{c}\cos \hat{\lambda }=k_1/\sqrt{{k_1}^2+{k_2}^2}\\ \sin \hat{\lambda }=k_2/\sqrt{{k_1}^2+{k_2}^2}\end{array}\right)& \left(\begin{array}{c}\sin \hat{\phi }=k_3\\ \cos \hat{\phi }=\sqrt{{k_1}^2+{k_2}^2}\end{array}\right)\end{array}---\left(32\right)$

由和确定λ和φ的表达式为

$\begin{array}{cc}\left(\begin{array}{c}\cos \lambda ={\tilde{k}}_1/\sqrt{{{\tilde{k}}_1}^2+{{\tilde{k}}_2}^2}\\ \sin \lambda ={\tilde{k}}_2/\sqrt{{{\tilde{k}}_1}^2+{{\tilde{k}}_2}^2}\end{array}\right)& \left(\begin{array}{c}\sin \phi ={\tilde{k}}_3\\ \cos \phi =\sqrt{{{\tilde{k}}_1}^2+{{\tilde{k}}_2}^2}\end{array}\right)\end{array}---\left(33\right)$

和σ的关系为，

$\hat{\sigma }=\sigma +\eta ---\left(34\right)$

其中，

$\left(\begin{array}{c}\sin \eta =\frac{\sin (\lambda -{\lambda }_P){\mathrm{cos\phi }}_P}{\cos \hat{\phi }}\\ \cos \eta =-\cos (A_P-\hat{\lambda })cos(\lambda -{\lambda }_P)+\sin (A_P-\hat{\lambda })\sin (\lambda -{\lambda }_P){\mathrm{sin\phi }}_P\end{array}\right)---\left(35\right)$

第四步，摄动方程建立

以飞行状态偏差量为状态变量，基于摄动理论建立换极坐标系中的运动摄动方程，

$\left(\begin{array}{c}\delta \stackrel{·}{\hat{r}}\\ \delta \stackrel{·}{\hat{\lambda }}\\ \delta \stackrel{·}{\hat{\phi }}\\ \delta \stackrel{·}{\hat{\sigma }}\\ \delta \stackrel{·}{\hat{\theta }}\\ \delta \stackrel{·}{\hat{V}}\end{array}\right)=\left(\begin{array}{cccccc}{A}_{11}& A_{12}& A_{13}& A_{14}& A_{15}& A_{16}\\ A_{21}& A_{22}& A_{23}& A_{24}& A_{25}& A_{26}\\ A_{31}& A_{32}& A_{33}& A_{34}& A_{35}& A_{36}\\ A_{41}& A_{42}& A_{43}& A_{44}& A_{45}& A_{46}\\ A_{51}& A_{52}& A_{53}& A_{54}& A_{55}& A_{56}\\ A_{61}& A_{62}& A_{63}& A_{64}& A_{65}& A_{66}\end{array}\right)\left(\begin{array}{c}\delta \hat{r}\\ \delta \hat{\lambda }\\ \delta \hat{\phi }\\ \delta \hat{\sigma }\\ \delta \hat{\theta }\\ \delta \hat{V}\end{array}\right)+\left(\begin{array}{c}{u}_1\\ u_2\\ u_3\\ u_4\\ u_5\\ u_6\end{array}\right)---\left(36\right)$

基于如下假设对A矩阵进行简化：

①不考虑地球自转；

②不考虑J2项；

③引入平衡滑翔假设，认为

④换极后弹道纵平面始终位于换极赤道附近，认为

进一步约去高阶小量，简化后A矩阵中不为零的元素为，

$\{\begin{array}{c}\begin{array}{ccccc}{A}_{15}=\hat{V}& A_{24}=\frac{\hat{V}\cos \hat{\sigma }}{\hat{r}}& A_{26}=\frac{\sin \hat{\sigma }}{\hat{r}}& A_{34}=-\frac{\hat{V}\sin \hat{\sigma }}{\hat{r}}& A_{36}=\frac{\cos \hat{\sigma }}{\hat{r}}\end{array}\\ \begin{array}{cccc}{A}_{43}=\frac{\hat{V}\sin \hat{\sigma }}{\hat{r}}& A_{46}=\frac{L\sin \hat{\upsilon }}{{\hat{V}}^2}& A_{56}=\frac{1}{\hat{r}}+\frac{\mu }{{\hat{V}}^2{\hat{r}}^2}& A_{65}=-\frac{\mu }{{\hat{r}}^2}\end{array}\end{array}---\left(37\right)$

第五步，摄动方程状态转移矩阵Φ(t_k,t)求解

先通过式(38)求解Φ(t_k,t)的伴随矩阵G(t,t_k)，

$\{\begin{array}{c}\frac{dG(t,t_k)}{dt}=-A^T\left(t\right)G(t,t_k)\\ G(t_k,t_k)=I_6\end{array}---\left(38\right)$

令

$\{\begin{array}{c}a=\frac{\hat{V}\sin \hat{\sigma }}{\hat{r}}\\ b=\sqrt{\frac{\mu }{{\hat{r}}^2}(\frac{1}{\hat{r}}+\frac{\mu }{{\hat{V}}^2{\hat{r}}^2})}=\frac{1}{\hat{V}{\hat{r}}^2}\sqrt{\mu ({\hat{V}}^2\hat{r}+\mu )}\end{array}---\left(39\right)$

解式(38)可得，

G₁₁＝1,G₁₂＝G₁₃＝G₁₄＝G₁₅＝G₁₆＝0(40)

G₂₂＝1,G₂₁＝G₂₃＝G₂₄＝G₂₅＝G₂₆＝0(41)

$\{\begin{array}{c}{G}_{32}=\frac{\cos \hat{\sigma }}{\sin \hat{\sigma }}-\frac{\cos \hat{\sigma }}{\sin \hat{\sigma }}\cos \left(a\tau \right)\\ G_{33}=\cos \left(a\tau \right)\\ G_{31}=G_{34}=G_{35}=G_{36}=0\end{array}---\left(42\right)$

$\{\begin{array}{c}{G}_{42}=-\frac{\cos \hat{\sigma }}{\sin \hat{\sigma }}\sin \left(a\tau \right)\\ G_{43}=\sin \left(a\tau \right)\\ G_{41}=G_{45}=G_{46}=0\\ G_{44}=1\end{array}---\left(43\right)$

$\{\begin{array}{c}{G}_{51}=-\frac{\hat{V}}{b}\sin \left(b\tau \right)\\ G_{52}=-\frac{\mu }{b^2{\hat{r}}^3\sin \hat{\sigma }}+(\frac{{\mathrm{\mu cos}}^2\hat{\sigma }}{b^2{\hat{r}}^3\sin \sigma }+\frac{\mu {\hat{V}}^{2}sin\hat{\sigma }{\cos }^2\hat{\sigma }}{b^4{\hat{r}}^5})\cos \left(a\tau \right)\\ -(\frac{\mu L\sin \hat{\upsilon }\cos \hat{\sigma }}{b^2{\hat{V}}^2{\hat{r}}^2\sin \hat{\sigma }}+\frac{\mu L\sin \hat{\sigma }\cos \hat{\sigma }\sin \hat{\upsilon }}{b^4{\hat{r}}^4})\sin \left(a\tau \right)+\frac{\mu \hat{V}{\cos }^2\hat{\sigma }}{b^3{\hat{r}}^4}\sin \left(a\tau \right)\sin \left(b\tau \right)\\ +\frac{\mu L\cos \hat{\sigma }\sin \hat{\upsilon }}{b^3\hat{V}{\hat{r}}^3}\cos \left(a\tau \right)\sin \left(b\tau \right)-\frac{\mu {\hat{V}}^2\sin \hat{\sigma }{\cos }^2\hat{\sigma }}{b^4{\hat{r}}^5}\cos \left(a\tau \right)\cos \left(b\tau \right)\\ +\frac{\mu L\sin \hat{\sigma }\cos \hat{\sigma }\sin \hat{\upsilon }}{b^4{\hat{r}}^4}\sin \left(a\tau \right)\cos \left(b\tau \right)\\ G_{53}=-(\frac{\mu \cos \hat{\sigma }}{b^2{\hat{r}}^3}+\frac{\mu {\hat{V}}^2{\sin }^2\hat{\sigma }\cos \hat{\sigma }}{b^4{\hat{r}}^5})\cos \left(a\tau \right)+(\frac{\mu L\sin \hat{\upsilon }}{b^2{\hat{V}}^2{\hat{r}}^2}+\frac{{\mathrm{\mu Lsin}}^2\hat{\sigma }\sin \hat{\upsilon }}{b^4{\hat{r}}^4})\sin \left(a\tau \right)\\ -\frac{\mu \hat{V}\sin \hat{\sigma }\cos \hat{\sigma }}{b^3{\hat{r}}^4}\sin \left(a\tau \right)\sin \left(b\tau \right)-\frac{\mu L\sin \hat{\sigma }\sin \hat{\upsilon }}{b^3{\hat{r}}^3\hat{V}}\cos \left(a\tau \right)\sin \left(b\tau \right)\\ +\frac{\mu {\hat{V}}^2{\sin }^2\hat{\sigma }\cos \hat{\sigma }}{b^4{\hat{r}}^5}\cos \left(a\tau \right)\cos \left(b\tau \right)-\frac{{\mathrm{\mu Lsin}}^2\hat{\sigma }\sin \hat{\upsilon }}{b^4{\hat{r}}^4}\sin \left(a\tau \right)\cos \left(b\tau \right)\\ G_{54}=\frac{\mu L\sin \hat{\upsilon }}{b^2{\hat{V}}^2{\hat{r}}^2}\\ G_{55}=0\\ G_{56}=\frac{\mu }{b{\hat{r}}^2}\sin \left(b\tau \right)\end{array}---\left(44\right)$

$\{\begin{array}{c}{G}_{61}=\frac{\hat{V}{\hat{r}}^2}{\mu }-\frac{\hat{V}{\hat{r}}^2}{\mu }\cos \left(b\tau \right)\\ G_{62}=\left(\cos \right(b\tau )-1)\left(\frac{\hat{V}{\cos }^2\hat{\sigma }}{b^2{\hat{r}}^2}\sin \right(a\tau )+\frac{L\cos \hat{\sigma }\sin \hat{\upsilon }}{b^2\hat{V}\hat{r}}\cos (a\tau \left)\right)\\ G_{63}=(1-\cos (b\tau \left)\right)\left(\frac{\hat{V}\sin \hat{\sigma }\cos \hat{\sigma }}{b^2{\hat{r}}^2}\sin \right(a\tau )+\frac{L\sin \hat{\sigma }\sin \hat{\upsilon }}{b^2\hat{V}\hat{r}}\cos (a\tau \left)\right)\\ G_{64}=G_{65}=0\\ G_{66}=\cos \left(b\tau \right)\end{array}---\left(45\right)$

再通过式(46)求解Φ(t_k,t)，

Φ(t_k,t)＝G^T(t,t_k)(46)

可得，

$\Phi (t_k,t)=\left(\begin{array}{cccccc}1& 0& 0& 0& G_{51}& G_{61}\\ 0& 1& G_{32}& G_{42}& G_{52}& G_{62}\\ 0& 0& G_{33}& G_{43}& G_{53}& G_{63}\\ 0& 0& 0& 1& G_{54}& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& G_{56}& G_{66}\end{array}\right)---\left(47\right)$

第六步，误差传播模型建立

一般坐标系中误差传播模型建立方法如下：

①根据式(48)建立换极坐标系中误差传播模型，

$\left(\begin{array}{c}\delta \hat{r}\\ \delta \hat{\lambda }\\ \delta \hat{\phi }\\ \delta \hat{\sigma }\\ \delta \hat{\theta }\\ \delta \hat{V}\end{array}\right)=\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right)\left(\begin{array}{c}\delta {\hat{r}}_0\\ \delta {\hat{\lambda }}_0\\ \delta {\hat{\phi }}_0\\ \delta {\hat{\sigma }}_0\\ \delta {\hat{\theta }}_0\\ \delta {\hat{V}}_0\end{array}\right)+{\int }_{t_0}^{t_k}\left(\begin{array}{cccccc}{\Phi }_{11}& {\Phi }_{12}& {\Phi }_{13}& {\Phi }_{14}& {\Phi }_{15}& {\Phi }_{16}\\ {\Phi }_{21}& {\Phi }_{22}& {\Phi }_{23}& {\Phi }_{24}& {\Phi }_{25}& {\Phi }_{26}\\ {\Phi }_{31}& {\Phi }_{32}& {\Phi }_{33}& {\Phi }_{34}& {\Phi }_{35}& {\Phi }_{36}\\ {\Phi }_{41}& {\Phi }_{42}& {\Phi }_{43}& {\Phi }_{44}& {\Phi }_{45}& {\Phi }_{46}\\ {\Phi }_{51}& {\Phi }_{52}& {\Phi }_{53}& {\Phi }_{54}& {\Phi }_{55}& {\Phi }_{56}\\ {\Phi }_{61}& {\Phi }_{62}& {\Phi }_{63}& {\Phi }_{64}& {\Phi }_{65}& {\Phi }_{66}\end{array}\right)\left(\begin{array}{c}{u}_1\\ u_2\\ u_3\\ u_4\\ u_5\\ u_6\end{array}\right)d\tau ---\left(48\right)$

②根据第三步中的坐标转换关系获得一般坐标系中误差传播模型。

至此，经过上述六步可最终建立一种基于摄动理论的滑翔弹道误差传播分析方法，该 方法具有“计算速度快、适应范围广、方法误差小”的特征，能满足滑翔弹道摄动因素误差 影响分析的精度要求。

本发明中，摄动因素包括飞行器本体建模误差和地球物理环境因素建模误差。在一种 具体实施方式中，摄动因素可为所述飞行器本体建模误差和地球物理环境因素建模误差中 的一种或多种。本领域技术人员能理解的，所述扰动引力是指地球对滑翔飞行器产生的扰 动引力。

与现有研究基础相比，本发明提出的方法具有以下优点：

1)首次提出一种针对高超声速滑翔飞行器滑翔弹道误差传播机理的分析方法，方法有 助于获得并深化对误差传播特性物理本质的机理认识，具有普遍适应性和广阔的工 程应用前景，为飞行器本体模型的精细化设计、弹道规划和制导方法误差的快速修 正、地球物理环境因素的影响补偿等提供理论基础和方法支撑。

2)首次将摄动理论应用于高超声速滑翔飞行器的弹道误差分析，以飞行状态偏差量为 状态变量建立运动摄动方程，将摄动因素处理为由运动摄动方程描述的动力系统的 外部激励。基于对该摄动问题属于正则摄动问题的判断，认为该动力系统的结构和 性质不受外部激励的影响，从而实现摄动因素和摄动方程动力系统的解耦。

3)首次将坐标变换的思路引入摄动方程状态转移矩阵的简化求解问题，通过建立一种 换极坐标系并推导换极坐标系中的飞行器动力学方程，使滑翔弹道再入平面始终处 于换极赤道面附近，从而可认为换极后的地心纬度近似为零，极大地简化了原问题 的求解。

4)通过求解伴随矩阵获得状态转移矩阵的解析解，将不同时刻齐次方程多次积分问题 简化为终端时刻单次逆向积分问题，使计算量大为减小。

5)方法的推导立足于飞行器动力学模型和弹道特性，可广泛适用于所有以相同动力学 模型描述的飞行器，不受飞行器特征参数改变的影响。方法所针对的摄动因素具有 普适性的数学描述且基于相同的处理手段，涵盖飞行器本体建模误差和地球物理环 境摄动因素，可适用于任意单一摄动因素或复杂摄动因素的影响分析。方法所建立 的误差传播模型具有解耦的零输入响应和零状态响应，分别对应飞行初态误差和飞 行过程中的摄动误差，可独立分析其影响。

6)该方法具有计算速度快、适用范围广、方法误差小等特征，可将误差影响的分析时 间缩短为弹道求差法的1/4，将方法误差在滑翔弹道终端飞行状态偏差量中的占比 控制在10％以内，将由方法误差引起的滑翔弹道终端位置偏差控制在2km以内； 方法适用于任意滑翔弹道，不出现奇点。

附图说明

图1为本发明换极坐标系及飞行状态量定义示意图；

图2为本发明滑翔弹道起点I和换极后极点P的位置关系示意图；

图3为本发明换极坐标系航迹偏航角的位置关系示意图；

图4为本发明换极坐标系中滑翔弹道示意图；

图5为本发明一般坐标系中滑翔弹道示意图；

图6为本发明沿滑翔弹道的扰动引力示意图；

图7为本发明飞行状态量及飞行控制量示意图；

图8为本发明单位扰动引力影响权系数示意图；

图9为本发明由误差传播模型计算得到的典型弹道终端位置偏差示意图；

图10为本发明由误差传播模型和弹道求差法计算得到的不同滑翔弹道终端位置偏差 示意图。

具体实施方式

下面结合具体实例，对本发明作进一步的说明：

基于CAV-H气动模型，以滑翔飞行弹道为研究对象，以飞行过程中扰动引力为单一摄 动因素进行仿真。飞行过程中弹道上各点的扰动引力采用1080阶球谐函数法进行计算。

仿真初始条件设置为：

(1)针对某典型滑翔弹道进行误差传播机理分析：①初始点状态参数：速度 V₀＝6500m/s，速度倾角θ₀＝0°，高度H₀＝80km，经度λ₀＝100°，地心纬度φ₀＝20°；② 结束点状态参数：速度V_k＝2000m/s，速度倾角θ_k＝0°，高度H_k＝25km，经度λ_k＝165°， 地心纬度φ_k＝50°；③终端结束时距目标点的待飞航程：S_togo＝100km。(2)针对具有不同 射向的系列滑翔弹道进行方法适应性分析：①初始点状态参数：速度V₀＝6500m/s，速度 倾角θ₀＝0°，高度H₀＝80km，经度λ₀＝100°，地心纬度φ₀＝0°；②结束点状态参数：速 度V_k＝2000m/s，速度倾角θ_k＝0°，高度H_k＝25km，终端经纬度见表1；③终端结束时距 目标点的待飞航程：S_togo＝100km。

表1不同射向滑翔弹道初始方位角及终点经纬度

序号 1 2 3 4 5 6 7 经度(deg) 100 138.93 149.25 154.74 157.85 159.49 160 地心纬度(deg) -60 -50 -40 -30 -20 -10 0 起始点方位角(deg) 180.01 152.21 137.93 125.27 113.27 101.57 90 序号 8 9 10 11 12 13 - 经度(deg) 100 138.93 149.25 154.74 157.85 159.49 - 地心纬度(deg) 60 50 40 30 20 10 - 起始点方位角(deg) 78.44 66.74 54.74 42.08 27.81 0 -

仿真计算机配置为：Intel(R)Core(TM)i5-3470CPU@3.20GHz，内存为3.46GB。软 件环境为WindowXP操作系统，计算程序基于VC++6.0开发。

其具体步骤如下：

第一步，换极坐标系建立

按如下方式建立换极坐标系：

①定义一个再入大圆弧平面作为换极赤道平面：1)对目标点确定的情况，将滑翔起 点和目标点地心矢径构成的再入大圆弧平面作为换极赤道平面；2)对于目标点未确定的 情况，根据滑翔起点位置及方位角确定的再入大圆弧平面作为换极赤道面。

②基于换极赤道平面定义换极坐标系O_E为地心，轴沿滑翔起点地心矢 径方向，轴在换极赤道面内垂直于轴指向目标点方向，轴与轴，轴构成右手系。 见图1。

第二步，换极坐标系中飞行器动力学模型建立

以高超声速滑翔飞行器滑翔弹道为研究对象，确定其飞行状态量为经度地心纬度航迹偏航角速度速度倾角和地心距建立以时间为自变量的动力学模型，

$\left(\begin{array}{c}\frac{d\hat{r}}{dt}=\hat{V}\sin \hat{\theta }\\ \frac{d\hat{\lambda }}{dt}=\frac{\hat{V}\cos \hat{\theta }\sin \hat{\sigma }}{\hat{r}\cos \hat{\phi }}\\ \frac{d\hat{\phi }}{dt}=\frac{\hat{V}\cos \hat{\theta }\cos \hat{\sigma }}{\hat{r}}\\ \frac{d\hat{\sigma }}{dt}=\frac{L\sin \upsilon }{\hat{V}\cos \hat{\theta }}+\frac{\hat{V}}{\hat{r}}\cos \hat{\theta }\tan \hat{\phi }\sin \hat{\sigma }-{\hat{g}}_{{\omega }_e}\frac{\cos \hat{\phi }\sin \hat{\sigma }}{\hat{V}\cos \hat{\theta }}+C_{\sigma }+{\tilde{C}}_{\sigma }\\ \frac{d\hat{\theta }}{dt}=\frac{L}{\hat{V}}\cos \upsilon +\frac{\hat{V}}{\hat{r}}\cos \hat{\theta }+\frac{{\hat{g}}_r\cos \hat{\theta }}{\hat{V}}+\frac{{\hat{g}}_{{\omega }_e}}{\hat{V}}(\cos \hat{\theta }\sin \hat{\phi }-\cos \hat{\sigma }\sin \hat{\theta }\cos \hat{\phi })+C_{\theta }+{\tilde{C}}_{\theta }\\ \frac{d\hat{V}}{dt}=-D+{\hat{g}}_r\sin \hat{\theta }+{\hat{g}}_{{\omega }_e}(\cos \hat{\sigma }\cos \hat{\theta }\cos \hat{\phi }+\sin \hat{\theta }\sin \hat{\phi })+{\tilde{C}}_V\end{array}\right)---\left(49\right)$

其中，C_σ、C_θ为哥氏加速度项，和为牵连加速度项，其表达式如下，

$\left(\begin{array}{c}{C}_{\sigma }=(2{\omega }_{ex}-2\tan \hat{\theta }({\omega }_{ey}\sin \hat{\sigma }+{\omega }_{ez}\cos \hat{\sigma }\left)\right)\\ {\tilde{C}}_{\sigma }=-\frac{\hat{r}}{\hat{V}\cos \hat{\theta }}({\omega }_{ex}{\omega }_{ey}\cos \hat{\sigma }-{\omega }_{ex}{\omega }_{ez}\sin \hat{\sigma })\\ C_{\theta }=2({\omega }_{ez}\sin \hat{\sigma }-{\omega }_{ey}\cos \hat{\sigma })\\ {\tilde{C}}_{\theta }=\frac{\hat{r}}{\hat{V}}[{\omega }_{ex}{\omega }_{ey}\sin \hat{\theta }\sin \hat{\sigma }+{\omega }_{ex}{\omega }_{ez}\sin \hat{\theta }\cos \hat{\sigma }+({\omega }_{ey}^2+{\omega }_{ez}^2)\cos \hat{\theta }]\\ {\tilde{C}}_V=\hat{r}[-{\omega }_{ex}{\omega }_{ey}\cos \hat{\theta }\sin \hat{\sigma }-{\omega }_{ex}{\omega }_{ez}\cos \hat{\theta }\cos \hat{\sigma }+({\omega }_{ey}^2+{\omega }_{ez}^2)\sin \hat{\theta }]\end{array}\right)---\left(50\right)$

其中，

$\left(\begin{array}{c}{\omega }_{ex}={\omega }_e(\cos \hat{\lambda }\cos \hat{\phi }{\mathrm{cos\phi }}_p{\mathrm{cosA}}_p+\sin \hat{\lambda }\cos \hat{\phi }{\mathrm{cos\phi }}_p{\mathrm{sinA}}_p+\sin \hat{\phi }{\mathrm{sin\phi }}_p)\\ {\omega }_{ex}={\omega }_e(-\sin \hat{\lambda }{\mathrm{cos\phi }}_p{\mathrm{cosA}}_p+\cos \hat{\lambda }{\mathrm{cos\phi }}_p{\mathrm{sinA}}_p)\\ {\omega }_{ez}={\omega }_e(-\cos \hat{\lambda }\sin \hat{\phi }{\mathrm{cos\phi }}_p{\mathrm{cosA}}_p-\sin \hat{\lambda }\sin \hat{\phi }{\mathrm{cos\phi }}_p{\mathrm{sinA}}_p+\cos \hat{\phi }{\mathrm{sin\phi }}_p)\end{array}\right)---\left(51\right)$

其中，ω_e为地球旋转加速度矢量，λ_p为换极后极点P的经度，φ_p为P的地心纬度，A_P为P的方位角，见图2。

第三步，一般坐标系与换极坐标系中飞行状态参数转换

根据换极坐标系定义，一般坐标系与换极坐标系中地心距、当地速度倾角及速度的定 义一致，

$r=\hat{r},\theta =\hat{\theta },V=\hat{V}---\left(52\right)$

定义

$\begin{array}{c}\left(\begin{array}{ccc}{\mathrm{cos\phi }}_0{\mathrm{cos\lambda }}_0& {\mathrm{cos\phi }}_0{\mathrm{sin\lambda }}_0& {\mathrm{sin\phi }}_0\\ -{\mathrm{sin\psi }}_0{\mathrm{sin\lambda }}_0-{\mathrm{cos\psi }}_0{\mathrm{sin\phi }}_0{\mathrm{cos\lambda }}_0& {\mathrm{sin\psi }}_0{\mathrm{cos\lambda }}_0-{\mathrm{cos\psi }}_0{\mathrm{sin\phi }}_0{\mathrm{sin\lambda }}_0& {\mathrm{cos\psi }}_0{\mathrm{cos\phi }}_0\\ {\mathrm{cos\psi }}_0{\mathrm{sin\lambda }}_0-{\mathrm{sin\psi }}_0{\mathrm{sin\phi }}_0{\mathrm{cos\lambda }}_0& -{\mathrm{cos\psi }}_0{\mathrm{cos\lambda }}_0-{\mathrm{sin\psi }}_0{\mathrm{sin\phi }}_0{\mathrm{sin\lambda }}_0& {\mathrm{sin\psi }}_0{\mathrm{cos\phi }}_0\end{array}\right)\\ \stackrel{\Delta }{=}\left(\begin{array}{ccc}{G}_{11}& G_{12}& G_{13}\\ G_{21}& G_{22}& G_{23}\\ G_{31}& G_{32}& G_{33}\end{array}\right).\end{array}---\left(53\right)$

$\{\begin{array}{c}{G}_{11}\cos \phi \cos \lambda +G_{12}\cos \phi \sin \lambda +G_{13}\sin \phi \stackrel{\Delta }{=}{k}_1\\ G_{21}\cos \phi \cos \lambda +G_{22}\cos \phi \sin \lambda +G_{23}\sin \phi \stackrel{\Delta }{=}{k}_2\\ G_{31}\cos \phi \cos \lambda +G_{32}\cos \phi \sin \lambda +G_{33}\sin \phi \stackrel{\Delta }{=}{k}_3\end{array}---\left(54\right)$

$\{\begin{array}{c}{G}_{11}\cos \hat{\phi }\cos \hat{\lambda }+G_{21}\cos \hat{\phi }\sin \hat{\lambda }+G_{31}\sin \hat{\phi }\stackrel{\Delta }{=}{\tilde{k}}_1\\ G_{12}\cos \hat{\phi }\cos \hat{\lambda }+G_{22}\cos \hat{\phi }\sin \hat{\lambda }+G_{32}\sin \hat{\phi }\stackrel{\Delta }{=}{\tilde{k}}_2\\ G_{13}\cos \hat{\phi }\cos \hat{\lambda }+G_{23}\cos \hat{\phi }\sin \hat{\lambda }+G_{33}\sin \hat{\phi }\stackrel{\Delta }{=}{\tilde{k}}_3\end{array}---\left(55\right)$

由λ和φ确定和的表达式为，

$\begin{array}{cc}\left(\begin{array}{c}\cos \hat{\lambda }=k_1/\sqrt{{k_1}^2+{k_2}^2}\\ \sin \hat{\lambda }=k_2/\sqrt{{k_1}^2+{k_2}^2}\end{array}\right)& \left(\begin{array}{c}\sin \hat{\phi }=k_3\\ \cos \hat{\phi }=\sqrt{{k_1}^2+{k_2}^2}\end{array}\right)\end{array}---\left(56\right)$

由和确定λ和φ的表达式为

$\begin{array}{cc}\left(\begin{array}{c}\cos \lambda ={\tilde{k}}_1/\sqrt{{{\tilde{k}}_1}^2+{{\tilde{k}}_2}^2}\\ \sin \lambda ={\tilde{k}}_2/\sqrt{{{\tilde{k}}_1}^2+{{\tilde{k}}_2}^2}\end{array}\right)& \left(\begin{array}{c}\sin \phi ={\tilde{k}}_3\\ \cos \phi =\sqrt{{{\tilde{k}}_1}^2+{{\tilde{k}}_2}^2}\end{array}\right)\end{array}---\left(57\right)$

确定σ的位置关系见图3，其表达式为，

$\hat{\sigma }=\sigma +\eta ---\left(58\right)$

其中，

$\left(\begin{array}{c}\sin \eta =\frac{\sin (\lambda -{\lambda }_P){\mathrm{cos\phi }}_P}{\cos \hat{\phi }}\\ \cos \eta =-\cos (A_P-\hat{\lambda })\cos (\lambda -{\lambda }_P)+\sin (A_P-\hat{\lambda })\sin (\lambda -{\lambda }_P){\mathrm{sin\phi }}_P\end{array}\right)---\left(59\right)$

第四步，摄动方程建立

以飞行状态偏差量为状态变量，基于摄动理论建立换极坐标系中的运动摄动方程，

$\left(\begin{array}{c}\delta \stackrel{·}{\hat{r}}\\ \delta \stackrel{·}{\hat{\lambda }}\\ \delta \stackrel{·}{\hat{\phi }}\\ \delta \stackrel{·}{\hat{\sigma }}\\ \delta \stackrel{·}{\hat{\theta }}\\ \delta \stackrel{·}{\hat{V}}\end{array}\right)=\left(\begin{array}{cccccc}{A}_{11}& A_{12}& A_{13}& A_{14}& A_{15}& A_{16}\\ A_{21}& A_{22}& A_{23}& A_{24}& A_{25}& A_{26}\\ A_{31}& A_{32}& A_{33}& A_{34}& A_{35}& A_{36}\\ A_{41}& A_{42}& A_{43}& A_{44}& A_{45}& A_{46}\\ A_{51}& A_{52}& A_{53}& A_{54}& A_{55}& A_{56}\\ A_{61}& A_{62}& A_{63}& A_{64}& A_{65}& A_{66}\end{array}\right)\left(\begin{array}{c}\delta \hat{r}\\ \delta \hat{\lambda }\\ \delta \hat{\phi }\\ \delta \hat{\sigma }\\ \delta \hat{\theta }\\ \delta \hat{V}\end{array}\right)+\left(\begin{array}{c}{u}_1\\ u_2\\ u_3\\ u_4\\ u_5\\ u_6\end{array}\right)---\left(60\right)$

基于如下假设对A矩阵进行简化：

①不考虑地球自转；

②不考虑J2项；

③引入平衡滑翔假设，认为

④换极后弹道纵平面始终位于换极赤道附近，认为

进一步约去高阶小量，简化后A矩阵中不为零的元素为，

$\{\begin{array}{c}\begin{array}{ccccc}{A}_{15}=\hat{V}& A_{24}=\frac{\hat{V}\cos \hat{\sigma }}{\hat{r}}& A_{26}=\frac{\sin \hat{\sigma }}{\hat{r}}& A_{34}=-\frac{\hat{V}\sin \hat{\sigma }}{\hat{r}}& A_{36}=\frac{\cos \hat{\sigma }}{\hat{r}}\end{array}\\ \begin{array}{cccc}{A}_{43}=\frac{\hat{V}\sin \hat{\sigma }}{\hat{r}}& A_{46}=\frac{L\sin \hat{\upsilon }}{{\hat{V}}^2}& A_{56}=\frac{1}{\hat{r}}+\frac{\mu }{{\hat{V}}^2{\hat{r}}^2}& A_{65}=-\frac{\mu }{{\hat{r}}^2}\end{array}\end{array}---\left(61\right)$

第五步，摄动方程状态转移矩阵Φ(t_k,t)求解

先通过式(62)求解Φ(t_k,t)的伴随矩阵G(t,t_k)，

$\{\begin{array}{c}\frac{dG(t,t_k)}{dt}=-A^T\left(t\right)G(t,t_k)\\ G(t_k,t_k)=I_6\end{array}---\left(62\right)$

令

$\{\begin{array}{c}a=\frac{\hat{V}\sin \hat{\sigma }}{\hat{r}}\\ b=\sqrt{\frac{\mu }{{\hat{r}}^2}(\frac{1}{\hat{r}}+\frac{\mu }{{\hat{V}}^2{\hat{r}}^2})}=\frac{1}{\hat{V}{\hat{r}}^2}\sqrt{\mu ({\hat{V}}^2\hat{r}+\mu )}\end{array}---\left(63\right)$

解式(62)可得，

G₁₁＝1,G₁₂＝G₁₃＝G₁₄＝G₁₅＝G₁₆＝0(64)

G₂₂＝1,G₂₁＝G₂₃＝G₂₄＝G₂₅＝G₂₆＝0(65)

$\{\begin{array}{c}{G}_{32}=\frac{\cos \hat{\sigma }}{\sin \hat{\sigma }}-\frac{\cos \hat{\sigma }}{\sin \hat{\sigma }}\cos \left(a\tau \right)\\ G_{33}=\cos \left(a\tau \right)\\ G_{31}=G_{34}=G_{35}=G_{36}=0\end{array}---\left(66\right)$

$\{\begin{array}{c}{G}_{42}=-\frac{\cos \hat{\sigma }}{\sin \hat{\sigma }}\sin \left(a\tau \right)\\ G_{43}=\sin \left(a\tau \right)\\ G_{41}=G_{45}=G_{46}=0\\ G_{44}=1\end{array}---\left(67\right)$

$\{\begin{array}{c}{G}_{51}=-\frac{\hat{V}}{b}\sin \left(b\tau \right)\\ G_{52}=-\frac{\mu }{b^2{\hat{r}}^3\sin \hat{\sigma }}+(\frac{{\mathrm{\mu cos}}^2\hat{\sigma }}{b^2{\hat{r}}^3\sin \sigma }+\frac{\mu {\hat{V}}^{2}sin\hat{\sigma }{\cos }^2\hat{\sigma }}{b^4{\hat{r}}^5})\cos \left(a\tau \right)\\ -(\frac{\mu L\sin \hat{\upsilon }\cos \hat{\sigma }}{b^2{\hat{V}}^2{\hat{r}}^2\sin \hat{\sigma }}+\frac{\mu L\sin \hat{\sigma }\cos \hat{\sigma }\sin \hat{\upsilon }}{b^4{\hat{r}}^4})\sin \left(a\tau \right)+\frac{\mu \hat{V}{\cos }^2\hat{\sigma }}{b^3{\hat{r}}^4}\sin \left(a\tau \right)\sin \left(b\tau \right)\\ +\frac{\mu L\cos \hat{\sigma }\sin \hat{\upsilon }}{b^3\hat{V}{\hat{r}}^3}\cos \left(a\tau \right)\sin \left(b\tau \right)-\frac{\mu {\hat{V}}^2\sin \hat{\sigma }{\cos }^2\hat{\sigma }}{b^4{\hat{r}}^5}\cos \left(a\tau \right)\cos \left(b\tau \right)\\ +\frac{\mu L\sin \hat{\sigma }\cos \hat{\sigma }\sin \hat{\upsilon }}{b^4{\hat{r}}^4}\sin \left(a\tau \right)\cos \left(b\tau \right)\\ G_{53}=-(\frac{\mu \cos \hat{\sigma }}{b^2{\hat{r}}^3}+\frac{\mu {\hat{V}}^2{\sin }^2\hat{\sigma }\cos \hat{\sigma }}{b^4{\hat{r}}^5})\cos \left(a\tau \right)+(\frac{\mu L\sin \hat{\upsilon }}{b^2{\hat{V}}^2{\hat{r}}^2}+\frac{{\mathrm{\mu Lsin}}^2\hat{\sigma }\sin \hat{\upsilon }}{b^4{\hat{r}}^4})\sin \left(a\tau \right)\\ -\frac{\mu \hat{V}\sin \hat{\sigma }\cos \hat{\sigma }}{b^3{\hat{r}}^4}\sin \left(a\tau \right)\sin \left(b\tau \right)-\frac{\mu L\sin \hat{\sigma }\sin \hat{\upsilon }}{b^3{\hat{r}}^3\hat{V}}\cos \left(a\tau \right)\sin \left(b\tau \right)\\ +\frac{\mu {\hat{V}}^2{\sin }^2\hat{\sigma }\cos \hat{\sigma }}{b^4{\hat{r}}^5}\cos \left(a\tau \right)\cos \left(b\tau \right)-\frac{{\mathrm{\mu Lsin}}^2\hat{\sigma }\sin \hat{\upsilon }}{b^4{\hat{r}}^4}\sin \left(a\tau \right)\cos \left(b\tau \right)\\ G_{54}=\frac{\mu L\sin \hat{\upsilon }}{b^2{\hat{V}}^2{\hat{r}}^2}\\ G_{55}=0\\ G_{56}=\frac{\mu }{b{\hat{r}}^2}\sin \left(b\tau \right)\end{array}---\left(68\right)$

$\{\begin{array}{c}{G}_{61}=\frac{\hat{V}{\hat{r}}^2}{\mu }-\frac{\hat{V}{\hat{r}}^2}{\mu }\cos \left(b\tau \right)\\ G_{62}=\left(\cos \right(b\tau )-1)\left(\frac{\hat{V}{\cos }^2\hat{\sigma }}{b^2{\hat{r}}^2}\sin \right(a\tau )+\frac{L\cos \hat{\sigma }\sin \hat{\upsilon }}{b^2\hat{V}\hat{r}}\cos (a\tau \left)\right)\\ G_{63}=(1-\cos (b\tau \left)\right)\left(\frac{\hat{V}\sin \hat{\sigma }\cos \hat{\sigma }}{b^2{\hat{r}}^2}\sin \right(a\tau )+\frac{L\sin \hat{\sigma }\sin \hat{\upsilon }}{b^2\hat{V}\hat{r}}\cos (a\tau \left)\right)\\ G_{64}=G_{65}=0\\ G_{66}=\cos \left(b\tau \right)\end{array}---\left(69\right)$

再通过式(70)求解Φ(t_k,t)，

Φ(t_k,t)＝G^T(t,t_k)(70)

可得，

$\Phi (t_k,t)=\left(\begin{array}{cccccc}1& 0& 0& 0& G_{51}& G_{61}\\ 0& 1& G_{32}& G_{42}& G_{52}& G_{62}\\ 0& 0& G_{33}& G_{43}& G_{53}& G_{63}\\ 0& 0& 0& 1& G_{54}& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& G_{56}& G_{66}\end{array}\right)---\left(71\right)$

第六步，误差传播模型建立

一般坐标系中误差传播模型建立方法如下：

①根据式(72)建立换极坐标系中误差传播模型，

$\left(\begin{array}{c}\delta \hat{r}\\ \delta \hat{\lambda }\\ \delta \hat{\phi }\\ \delta \hat{\sigma }\\ \delta \hat{\theta }\\ \delta \hat{V}\end{array}\right)=\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right)\left(\begin{array}{c}\delta {\hat{r}}_0\\ \delta {\hat{\lambda }}_0\\ \delta {\hat{\phi }}_0\\ \delta {\hat{\sigma }}_0\\ \delta {\hat{\theta }}_0\\ \delta {\hat{V}}_0\end{array}\right)+{\int }_{t_0}^{t_k}\left(\begin{array}{cccccc}{\Phi }_{11}& {\Phi }_{12}& {\Phi }_{13}& {\Phi }_{14}& {\Phi }_{15}& {\Phi }_{16}\\ {\Phi }_{21}& {\Phi }_{22}& {\Phi }_{23}& {\Phi }_{24}& {\Phi }_{25}& {\Phi }_{26}\\ {\Phi }_{31}& {\Phi }_{32}& {\Phi }_{33}& {\Phi }_{34}& {\Phi }_{35}& {\Phi }_{36}\\ {\Phi }_{41}& {\Phi }_{42}& {\Phi }_{43}& {\Phi }_{44}& {\Phi }_{45}& {\Phi }_{46}\\ {\Phi }_{51}& {\Phi }_{52}& {\Phi }_{53}& {\Phi }_{54}& {\Phi }_{55}& {\Phi }_{56}\\ {\Phi }_{61}& {\Phi }_{62}& {\Phi }_{63}& {\Phi }_{64}& {\Phi }_{65}& {\Phi }_{66}\end{array}\right)\left(\begin{array}{c}{u}_1\\ u_2\\ u_3\\ u_4\\ u_5\\ u_6\end{array}\right)d\tau ---\left(72\right)$

②以飞行过程中扰动引力为单一摄动因素的换极坐标系误差传播模型为，

$\left(\begin{array}{c}\delta \hat{r}\\ \delta \hat{\lambda }\\ \delta \hat{\phi }\\ \delta \hat{\sigma }\\ \delta \hat{\theta }\\ \delta \hat{V}\end{array}\right)={\int }_{t_0}^{t_k}\left(\begin{array}{cccccc}{\Phi }_{11}& {\Phi }_{12}& {\Phi }_{13}& {\Phi }_{14}& {\Phi }_{15}& {\Phi }_{16}\\ {\Phi }_{21}& {\Phi }_{22}& {\Phi }_{23}& {\Phi }_{24}& {\Phi }_{25}& {\Phi }_{26}\\ {\Phi }_{31}& {\Phi }_{32}& {\Phi }_{33}& {\Phi }_{34}& {\Phi }_{35}& {\Phi }_{36}\\ {\Phi }_{41}& {\Phi }_{42}& {\Phi }_{43}& {\Phi }_{44}& {\Phi }_{45}& {\Phi }_{46}\\ {\Phi }_{51}& {\Phi }_{52}& {\Phi }_{53}& {\Phi }_{54}& {\Phi }_{55}& {\Phi }_{56}\\ {\Phi }_{61}& {\Phi }_{62}& {\Phi }_{63}& {\Phi }_{64}& {\Phi }_{65}& {\Phi }_{66}\end{array}\right)\left(\begin{array}{c}0\\ 0\\ 0\\ u_4\\ u_5\\ u_6\end{array}\right)d\tau ---\left(73\right)$

其中，摄动因素u₄、u₅和u₆为，

$\{\begin{array}{c}{u}_4=-\delta {\hat{g}}_{{\omega }_e}\frac{\cos \phi \sin \sigma }{V\cos \theta }\\ u_5=\frac{\delta {\hat{g}}_r\cos \theta }{V}+\frac{\delta {\hat{g}}_{{\omega }_e}}{V}(\cos \theta \sin \phi -\cos \sigma \sin \theta \cos \phi )\\ u_6=\delta {\hat{g}}_r\sin \theta +\delta {\hat{g}}_{{\omega }_e}(\cos \sigma \cos \theta \cos \phi +\sin \theta \sin \phi )\end{array}---\left(74\right)$

由扰动引力引起的位置偏差为，

$\left(\begin{array}{c}\delta r={\int }_0^{t_k}{C}_{\delta {\hat{g}}_{{\omega }_e}}^{\delta r}(t_k-\tau )\delta {\hat{g}}_{{\omega }_e}\left(\tau \right)d\tau +{\int }_0^{t_k}{C}_{\delta {\hat{g}}_r}^{\delta r}(t_k-\tau )\delta {\hat{g}}_r\left(\tau \right)d\tau \\ \delta \lambda ={\int }_0^{t_k}{C}_{\delta {\hat{g}}_{{\omega }_e}}^{\delta \lambda }(t_k-\tau )\delta {\hat{g}}_{{\omega }_e}\left(\tau \right)d\tau +{\int }_0^{t_k}{C}_{\delta {\hat{g}}_r}^{\delta \lambda }(t_k-\tau )\delta {\hat{g}}_r\left(\tau \right)d\tau \\ \delta \phi ={\int }_0^{t_k}{C}_{\delta {\hat{g}}_{{\omega }_e}}^{\delta \phi }(t_k-\tau )\delta {\hat{g}}_{{\omega }_e}\left(\tau \right)d\tau +{\int }_0^{t_k}{C}_{\delta {\hat{g}}_r}^{\delta \phi }(t_k-\tau )\delta {\hat{g}}_r\left(\tau \right)d\tau \end{array}\right)---\left(75\right)$

其中，$C_{\delta {\hat{g}}_{{\omega }_e}}^{\delta r}(t_k-\tau ),C_{\delta {\hat{g}}_r}^{\delta r}(t_k-\tau ),C_{\delta {\hat{g}}_{{\omega }_e}}^{\delta \lambda }(t_k-\tau ),C_{\delta {\hat{g}}_r}^{\delta \lambda }(t_k-\tau ),C_{\delta {\hat{g}}_{{\omega }_e}}^{\delta \phi }(t_k-\tau )$和$C_{\delta {\hat{g}}_r}^{\delta \phi }(t_k-\tau )$表征单位扰动引力对滑翔弹道终端位置影响的权系数，代入状态转移矩阵Φ(t_k,t)中各元素 可得其表达式为，

③根据第三步中的坐标转换关系获得一般坐标系中误差传播模型。

图4～图5给出了换极坐标系及一般坐标系中考虑扰动引力和不考虑扰动引力的滑翔弹 道曲线。计算结果表明，在一般坐标系中，飞行过程中扰动引力引起的终端高度偏差为 4.005m，经度偏差(东向偏差)为-0.2634°，地心纬度偏差(北向偏差)为-0.0627°；在换 极坐标系中，飞行过程中扰动引力引起的终端高度偏差为4.005m，经度偏差(东向偏差) 为-0.1773°，地心纬度偏差(北向偏差)为-0.0480°。可粗略估算，由滑翔段扰动引力引起 的总位置偏差可达几十公里。

图6给出了通过1080阶球谐函数法计算得到的沿典型滑翔弹道的扰动引力值。由图6 可见，由于滑翔弹道飞行高度低且覆盖区域广，沿飞行弹道的扰动引力变化较为复杂。其 中，沿地心矢径方向的扰动引力分量最大可达42.7mgal，最小为-85.5mgal；沿地球自转方 向的扰动引力分量最大可达85.1mgal，最小为5.4mgal。

图7所示为飞行状态量及飞行控制量随飞行时间的变化曲线。由图7可见，换极后的 地心纬度及速度倾角都在0°附近，表明第四步中的简化假设合理。

图8为单位扰动引力影响权系数示意图。由图8可见，扰动引力对终端位置的影响不 具有线性变化关系，经多次增大-减小-增大的变化以后，在弹道终端趋近于0。经粗略估算， 单位经度偏差或单位纬度偏差造成的横向距离偏差约110km左右。对比横侧向权系数和纵 向权系数可知，扰动引力对横侧向弹道终端位置的影响在几到几十千米，远大于对纵向的 影响。除受速度和地心距影响外，纵向位置偏差主要受速度方位角的影响，横侧向位置偏 差主要受倾侧角和速度方位角的影响。

图9为由误差传播模型计算得到的等时终端位置偏差示意图，表2为通过弹道求差法 计算得到的等时终端位置偏差。图9与表2计算结果较为接近。

表2扰动引力引起的终端位置偏差

图10给出了由误差传播模型和弹道求差法计算得到的不同射向滑翔弹道终端位置偏 差。以弹道求差法计算得到的两种弹道的等时偏差为标准，可计算本发明的方法误差。表 3给出了将误差传播模型应用于不同滑翔弹道终端位置偏差分析的方法误差。结果表明， 本发明的方法误差小于总偏差量的10％。

表3误差传播模型方法误差

初始方位角 180.01 152.21 137.93 125.27 113.27 101.57 90 高度误差(m) -11.73 7.01 -11.16 -8.23 1.73 -3.19 -2.35 经度(deg) 0.010 -0.010 -0.006 -0.006 -0.011 -0.018 -0.011 地心纬度(deg) 0.009 0.017 0.008 -0.094 0.004 0.003 0.013 初始方位角 78.44 66.74 54.74 42.08 27.81 0 - 高度误差(m) 2.99 -5.78 -5.01 0.37 -3.69 16.79 - 经度(deg) -0.016 -0.011 -0.013 -0.001 -0.019 -0.014 - 地心纬度(deg) -0.005 -0.032 -0.005 -0.002 -0.001 -0.009 -

表4为针对某典型滑翔弹道利用本发明和弹道求差法进行终端状态偏差分析所需要的 计算时间。由表4可见，误差传播分析方法的主要计算量为各点扰动引力计算，约耗时 9561.204秒，而误差传播模型计算仅需0.921秒；其计算总时间(滑翔弹道误差传播分析 总时间)约为弹道求差法的1/4。

表4误差传播方法和弹道求差法计算时间对比

综合上述仿真结果可获得以下结论：

1)由本发明确定的方法进行滑翔弹道终端状态偏差量分析，由方法误差引起的终端位 置偏差在公里量级，方法误差在偏差量中的占比可控制在10％以内，能够满足滑 翔弹道误差分析的精度要求；

2)由本发明确定的方法进行滑翔弹道终端状态偏差量分析，其计算时间仅为弹道求差 法的1/4，其中误差传播模块的计算时间小于1s；本发明以牺牲可容忍的小量精度 为代价提高计算速度，具有广阔的工程应用前景；

3)由本发明确定的误差传播模型，可半解析地反映终端状态偏差量与飞行状态量和飞 行控制量之间的数学关系，有助于加深机理认识和建立补偿机制；

4)本发明能够针对多种摄动因素、可适应各种条件下滑翔弹道的误差传播分析，具有 适应范围广的特征。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人 员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、 等同替换、改进等，均应包含在本发明的保护范围之内。