CN105184109A - 扰动引力作用下弹道助推段状态偏差解析计算方法 - Google Patents

Abstract

本发明提供了一种扰动引力作用下弹道助推段状态偏差解析计算方法，具体是以弹道导弹为研究对象，针对弹道导弹助推段状态偏差快速求解的问题，具体方法是：首先，根据发射任务设计并生成标准弹道；其次，运用摄动思想导出的状态偏差解析法求解助推段每一点的状态偏差；之后，基于助推段弹道每一点的状态偏差计算出对应点的视加速度偏差，并将该偏差视为扰动引力的高阶项，同时进行补偿；最终，运用牛顿迭代法对状态偏差结果进行修正。该方法能够实现沿任意飞行弹道助推段状态偏差的快速计算，其计算精度能够满足弹道计算的要求，计算速度远优于现有方法，为实现弹道快速机动发射奠定了基础。

Description

扰动引力作用下弹道助推段状态偏差解析计算方法

技术领域

本发明涉及飞行器动力学建模领域，具体涉及一种扰动引力作用下弹道助推段状态偏差解析计算方法。

背景技术

不断提升射前机动快速发射性能及命中精度是我国新一代弹道导弹发展的必然趋势。目前来看，制约我国弹道导弹命中精度的主要因素有制导工具误差和制导方法误差，随着惯性测量系统硬件水平的提高，可有效修正部分制导工具误差，使得弹道导弹制导工具误差逐渐降低。同时，制导方法误差的影响日益突出，而扰动引力是影响制导方法误差的主要因素。因此，建立快速、精确的扰动引力作用下弹道导弹助推段关机点状态偏差解析计算方法具有重大的军事意义和工程价值。

目前，导弹飞行力学领域的学者对大气层内导弹误差传播模型也有较深入的研究，但主要集中在对制导工具误差进行建模分析，而对引力模型、大气模型等模型误差对导弹状态偏差影响的研究则很少，通常只是采用弹道求差法分析扰动引力对助推段关机点状态偏差的影响特性。弹道求差法虽然能够精确的求得扰动引力作用下助推段关机点状态偏差量，但该方法基于弹道积分，计算耗时较长，需要的存储量大，不利于新形势下导弹发射快速性的要求。

因此，亟待建立一种面向快速机动发射应用的弹道导弹助推段状态偏差计算方法，该问题存在的难点为：一是在保证求解精度的前提下将复杂的导弹非线性动力学偏差方程合理简化为可导出完整解析解表达式的形式；二是需要考虑扰动引力与视加速度之间的耦合特性对助推段状态偏差的影响，扰动引力会导致导弹飞行状态的改变，导弹状态的改变又使得视加速度发生改变，因此，必须对视加速度偏差引起的导弹助推段状态偏差进行修正。

发明内容

本发明目的在于提供一种扰动引力作用下弹道助推段状态偏差解析计算方法，该方法首先需在飞行任务要求下设计一条标准弹道；其次需运用摄动思想导出的状态偏差解析法求解助推段每一点的状态偏差；再次需要基于助推段弹道每一点的状态偏差计算出对应点的视加速度偏差，并将该偏差理解为扰动引力的高阶项，同时进行补偿；最后运用牛顿迭代法对状态偏差进行迭代修正。具体技术方案如下：

一种扰动引力作用下弹道助推段状态偏差解析计算方法，包括以下步骤：第一步，助推段标准弹道设计；第二步，状态偏差模型建立；第三步，助推段扰动量拟合；第四步，视加速度偏差量求解；第五步，解析模型迭代修正。

以上技术方案中优选的，所述第一步中助推段标准弹道设计具体为：按照发射点、目标点、某一型号导弹的总体参数、大气模型以及地球引力模型的任务条件设计一条满足条件的弹道，并按照一定时间间隔保存助推段发射系中时间、速度、位置、质量、程序角五项状态值。

以上技术方案中优选的，所述第二步中状态偏差模型建立具体为：

标准弹道计算模型表示为表达式(1)，详情如下：

$\left\{\begin{array}{c}{\stackrel{\·}{\stackrel{\‾}{V}}}^{*}={\stackrel{\·}{\stackrel{\‾}{W}}}^{*}+g^{*}\left({\stackrel{\‾}{\mathrm{\ρ}}}^{*}\right)\\ {\stackrel{\·}{\stackrel{\‾}{\mathrm{\ρ}}}}^{*}={\stackrel{\‾}{V}}^{*}\end{array}\right.---\left(1\right);$

其中，为标准弹道视加速度，为标准弹道引力加速度；

若考虑扰动引力对导弹运动的影响，则实际飞行弹道与标准弹道的等时变分表示为表达式(2)，详情如下：

$\left\{\begin{array}{c}\mathrm{\δ}\stackrel{\·}{\stackrel{\‾}{V}}=\stackrel{\·}{\stackrel{\‾}{W}}-{\stackrel{\·}{\stackrel{\‾}{W}}}^{*}+\stackrel{\‾}{g}\left(\stackrel{\‾}{\mathrm{\ρ}}\right)-{\stackrel{\‾}{g}}^{*}\left({\stackrel{\‾}{\mathrm{\ρ}}}^{*}\right)\\ \mathrm{\δ}\stackrel{\·}{\stackrel{\‾}{\mathrm{\ρ}}}=\mathrm{\δ}\stackrel{\‾}{V}\end{array}\right.---\left(2\right);$

其中，为实际弹道视加速度，为实际弹道引力加速度；

将表达式(2)改写为矩阵形式，同时经过小偏差处理并略去高阶小量后可得表达式(3)，详情如下：

$\left(\begin{array}{c}\mathrm{\δ}\stackrel{\·}{\stackrel{\‾}{V}}\\ \mathrm{\δ}\stackrel{\·}{\stackrel{\‾}{\mathrm{\ρ}}}\end{array}\right)=T\left(\begin{array}{c}\mathrm{\δ}\stackrel{\‾}{V}\\ \mathrm{\δ}\stackrel{\‾}{\mathrm{\ρ}}\end{array}\right)+\left(\begin{array}{c}\mathrm{\δ}\stackrel{\‾}{F}\\ 0\end{array}\right)---\left(3\right);$

其中，为扰动量，且有为扰动引力，为视加速度偏差；式中T的表达式详见表达式(4)：

表达式(4)中A₁、A₂、A₃、A₄、A₅以及A₆六者各自的表达式详见(5)：

$\left\{\begin{array}{c}{A}_1=-n_b^2(1-3\frac{x^2}{r^2})\\ A_2=3n_b^2\frac{x}{r}\frac{R_0+y}{r}\\ A_3=3n_b^2\frac{x}{r}\·\frac{z}{r}\\ A_4=-n_b^2(1-3\frac{{(R_0+y)}^2}{r^2})\\ A_5=3n_b^2\frac{z}{r}\·\frac{R_0+y}{r}\\ A_6=-n_b^2(1-3\frac{z^2}{r^2})\end{array}\right.---\left(5\right);$

经过合理简化，采用伴随方程求得表达式(3)的状态转移矩阵解析解为表达式(7)，详情如下：

$\begin{array}{c}\mathrm{\Φ}(t_k,t)=\mathrm{\Φ}\left(\mathrm{\τ}\right)\\ =\left[\begin{array}{cccccc}{\mathrm{cosn}}_b\mathrm{\τ}& {\mathrm{\Φ}}_{12}\left(\mathrm{\τ}\right)& 0& -n_b\sin n_b\mathrm{\τ}& 0& 0\\ {\mathrm{\Φ}}_{21}\left(\mathrm{\τ}\right)& \cosh \sqrt{2}{n}_b\mathrm{\τ}& 0& 0& \sqrt{2}{n}_b\sinh \sqrt{2}{n}_b\mathrm{\τ}& 0\\ 0& 0& \cos n_b\mathrm{\τ}& 0& 0& -n_b\sin n_b\mathrm{\τ}\\ \frac{\sin n_b\mathrm{\τ}}{n_b}& 0& 0& {\mathrm{cosn}}_b\mathrm{\τ}& 0& 0\\ 0& \frac{\sinh \sqrt{2}{n}_b\mathrm{\τ}}{\sqrt{2}{n}_b}& 0& 0& \cosh \sqrt{2}{n}_b\mathrm{\τ}& 0\\ 0& 0& \frac{\sin n_b\mathrm{\τ}}{n_b}& 0& 0& \cos n_b\mathrm{\τ}\end{array}\right]\end{array}---\left(7\right);$

其中Φ₁₂(τ)详见表达式(8)：

${\mathrm{\Φ}}_{12}\left(\mathrm{\τ}\right)={\mathrm{\Φ}}_{21}\left(\mathrm{\τ}\right)=\frac{{\left(n_b\mathrm{\τ}\right)}^2}{3!}(1+\frac{2}{3}\frac{\mathrm{\τ}}{t_k}+\frac{1}{6}\frac{{\mathrm{\τ}}^2}{t_k^2})---\left(8\right);$

根据表达式(7)即可推得助推段扰动引力引起的状态偏差半解析表达式(9)，如下：

$\begin{array}{c}{\mathrm{\δv}}_{ax}\left(t_k\right)={\∫}_0^{t_k}{\mathrm{\Φ}}_{11}(t_k-\mathrm{\τ}){\mathrm{\δF}}_x\left(\mathrm{\τ}\right)d\mathrm{\τ}+{\∫}_0^{t_k}{\mathrm{\Φ}}_{12}(t_k-\mathrm{\τ}){\mathrm{\δF}}_y\left(\mathrm{\τ}\right)d\mathrm{\τ}\\ {\mathrm{\δv}}_{ay}\left(t_k\right)={\∫}_0^{t_k}{\mathrm{\Φ}}_{22}(t_k-\mathrm{\τ}){\mathrm{\δF}}_y\left(\mathrm{\τ}\right)d\mathrm{\τ}+{\∫}_0^{t_k}{\mathrm{\Φ}}_{21}(t_k-\mathrm{\τ}){\mathrm{\δF}}_x\left(\mathrm{\τ}\right)d\mathrm{\τ}\\ {\mathrm{\δv}}_{az}\left(t_k\right)={\∫}_0^{t_k}{\mathrm{\Φ}}_{33}(t_k-\mathrm{\τ}){\mathrm{\δF}}_z\left(\mathrm{\τ}\right)d\mathrm{\τ}\\ {\mathrm{\δx}}_a\left(t_k\right)={\∫}_0^{t_k}{\mathrm{\Φ}}_{41}(t_k-\mathrm{\τ}){\mathrm{\δF}}_x\left(\mathrm{\τ}\right)d\mathrm{\τ}\\ {\mathrm{\δy}}_a\left(t_k\right)={\∫}_0^{t_k}{\mathrm{\Φ}}_{52}(t_k-\mathrm{\τ}){\mathrm{\δF}}_y\left(\mathrm{\τ}\right)d\mathrm{\τ}\\ {\mathrm{\δz}}_a\left(t_k\right)={\∫}_0^{t_k}{\mathrm{\Φ}}_{63}(t_k-\mathrm{\τ}){\mathrm{\δF}}_z\left(\mathrm{\τ}\right)d\mathrm{\τ}\end{array}---\left(9\right).$

以上技术方案中优选的，所述合理简化的过程具体为：

由表达式(5)可知，A₁、A₂、A₃、A₄、A₅以及A₆六者均为小量，其中A₁、A₄以及A₆较大，其数量级与相当，而A₂比小一个数量级，A₃和A₅则比小两三个数量级，即有

$\left\{\begin{array}{c}\frac{R_0+y}{r}\&ap;1\\ \frac{z^2}{r^3}<<1\\ \frac{x^2}{r^2}<<1\end{array}\right.;$

因此，将表达式(5)改写为表达式(6)的形式，详情如下：

$\left\{\begin{array}{c}{A}_1=A_6=-n_b^2\\ A_4=2n_b^2\\ A_2=\frac{3x}{r}{n}_b^2\\ A_3=A_5=0\end{array}\right.---\left(6\right);$

同时，考虑到n_b在整个助推段过程中的变化范围小于1％，因此，积分时设n_b为常数。

以上技术方案中优选的，所述第三步中助推段扰动量拟合具体为：

采用最小二乘法将助推段扰动量拟合为关于时间的多项式函数，且发射惯性系中扰动量三分量分别进行拟合，具体是：设观测方程为表达式(10)：

Z_n×1＝H_n×(m+1)λ_(m+1)×1+ε(10)；

其中，Z_n×1为观测向量，n为拟合点的个数；H_n×(m+1)为系数矩阵，m代表了多项式拟合次数；λ_(m+1)×1为待定参数向量；ε为随机误差向量；

拟合多项式采用经验公式(11)则可得系数矩阵和扰动量拟合系数为表述式(12)和(13)，详情如下：

$\mathrm{\&delta;}F\left(t\right)=\underset{i=0}{\overset{n}{\&Sigma;}}{\mathrm{\&lambda;}}_i{\left(\frac{t}{10}\right)}^i---\left(11\right);$

λ＝(H^TH)^-1H^TZ(13)；

当拟合多项式取5阶时，可得扰动量拟合表达式(14)，如下：

$\left\{\begin{array}{c}{\mathrm{\&delta;F}}_x\left(t\right)={\mathrm{\&lambda;}}_{x0}+{\mathrm{\&lambda;}}_{x1}(t/10)+{\mathrm{\&lambda;}}_{x2}{(t/10)}^2+{\mathrm{\&lambda;}}_{x3}{(t/10)}^3+{\mathrm{\&lambda;}}_{x4}{(t/10)}^4+{\mathrm{\&lambda;}}_{x5}{(t/10)}^5\\ {\mathrm{\&delta;F}}_y\left(t\right)={\mathrm{\&lambda;}}_{y0}+{\mathrm{\&lambda;}}_{y1}(t/10)+{\mathrm{\&lambda;}}_{y2}{(t/10)}^2+{\mathrm{\&lambda;}}_{y3}{(t/10)}^3+{\mathrm{\&lambda;}}_{y4}{(t/10)}^4+{\mathrm{\&lambda;}}_{y5}{(t/10)}^5\\ {\mathrm{\&delta;F}}_z\left(t\right)={\mathrm{\&lambda;}}_{z0}+{\mathrm{\&lambda;}}_{z1}(t/10)+{\mathrm{\&lambda;}}_{z2}{(t/10)}^2+{\mathrm{\&lambda;}}_{z3}{(t/10)}^3+{\mathrm{\&lambda;}}_{z4}{(t/10)}^4+{\mathrm{\&lambda;}}_{z5}{(t/10)}^5\end{array}\right.---\left(14\right);$

将表达式(14)带入助推段扰动引力引起的状态偏差半解析表达式(9)导出助推段扰动引力引起的状态偏差完整解析表达式(15)、(16)、(17)、(18)、(19)以及(20)，详情如下：

$\begin{array}{c}{\mathrm{\&delta;v}}_{ax}\left(t_k\right)=\frac{1}{n_b}{\mathrm{\&lambda;}}_{x0}\sin \left(n_b{t}_k\right)+\frac{1}{{n_b}^2}{\mathrm{\&lambda;}}_{x1}(1-\cos (n_b{t}_k)+\frac{2}{{n_b}^3}{\mathrm{\&lambda;}}_{x2}(n_b{t}_k-\sin \left(n_b{t}_k\right))+\\ \frac{3}{{n_b}^4}{\mathrm{\&lambda;}}_{x3}(-2+n_b^2{t}_k^2+2\cos (n_b{t}_k\left)\right)+\frac{4}{{n_b}^5}{\mathrm{\&lambda;}}_{x4}(-6n_b{t}_k+n_b^3{t}_k^3+6\sin (n_b{t}_k\left)\right)+\\ \frac{5}{{n_b}^6}{\mathrm{\&lambda;}}_{x5}(24-12n_b^2{t}_k^2+n_b^4{t}_k^4-24\cos (n_b{t}_k\left)\right)+\frac{1}{90720}\&CenterDot;\{n_b^2{t}_k^3\&CenterDot;\&lsqb;8064{\mathrm{\&lambda;}}_{y0}+\\ t_k\&CenterDot;(1848{\mathrm{\&lambda;}}_{y1}+(696{\mathrm{\&lambda;}}_{y2}+333{\mathrm{\&lambda;}}_{y3}{t}_k+184333{\mathrm{\&lambda;}}_{y4}{t}_k^2+112{\mathrm{\&lambda;}}_{y5}{t}_k^3\left)\right)\&rsqb;\}\end{array}---\left(15\right);$

$\begin{array}{c}{\mathrm{\&delta;v}}_{ay}\left(t_k\right)=\frac{1}{\sqrt{2}{n}_b}{\mathrm{\&lambda;}}_{y0}\sinh \left(\sqrt{2}{n}_b{t}_k\right)+\frac{1}{2{n_b}^2}{\mathrm{\&lambda;}}_{y1}\left(\cosh \right(\sqrt{2}{n}_b{t}_k)-1)+\frac{2}{2{n_b}^3}{\mathrm{\&lambda;}}_{y2}(-2n_b{t}_k+\\ \sqrt{2}\sinh \left(\sqrt{2}{n}_b{t}_k\right))+\frac{3}{2{n_b}^4}{\mathrm{\&lambda;}}_{y3}(-1-n_b^2{t}_k^2+\cosh (\sqrt{2}{n}_b{t}_k\left)\right)+\frac{1}{{n_b}^5}{\mathrm{\&lambda;}}_{y4}(-2n_b{t}_k(3+\\ n_b^2{t}_k^2)+3\sqrt{2}\sinh \left(\right(\sqrt{2}{n}_b{t}_k\left)\right)-\frac{5}{{n_b}^6}{\mathrm{\&lambda;}}_{y5}(6+12n_b^2{t}_k^2+n_b^4{t}_k^4-6\cosh (\sqrt{2}{n}_b{t}_k\left)\right)+\\ \frac{1}{90720}\&CenterDot;\{n_b^2{t}_k^3\&CenterDot;\&lsqb;8064{\mathrm{\&lambda;}}_{x0}+t_k\&CenterDot;(1848{\mathrm{\&lambda;}}_{x1}+(696{\mathrm{\&lambda;}}_{x2}+333{\mathrm{\&lambda;}}_{x3}{t}_k+184333{\mathrm{\&lambda;}}_{x4}{t}_k^2+\\ 112{\mathrm{\&lambda;}}_{x5}{t}_k^3\left)\right)\&rsqb;\}\end{array}---\left(16\right);$

$\begin{array}{c}{\mathrm{\&delta;v}}_{az}\left(t_k\right)=\frac{1}{n_b}{\mathrm{\&lambda;}}_{z0}\sin \left(n_b{t}_k\right)+\frac{1}{{n_b}^2}{\mathrm{\&lambda;}}_{z1}(1-\cos (n_b{t}_k)+\frac{2}{{n_b}^3}{\mathrm{\&lambda;}}_{z2}(n_b{t}_k-\sin \left(n_b{t}_k\right))+\\ \frac{3}{{n_b}^4}{\mathrm{\&lambda;}}_{z3}(-2+n_b^2{t}_k^2+2\cos (n_b{t}_k\left)\right)+\frac{4}{{n_b}^5}{\mathrm{\&lambda;}}_{z4}(-6n_b{t}_k+n_b^3{t}_k^3+6\sin (n_b{t}_k\left)\right)+\\ \frac{5}{{n_b}^6}{\mathrm{\&lambda;}}_{z5}(24-12n_b^2{t}_k^2+n_b^4{t}_k^4-24\cos (n_b{t}_k\left)\right)\end{array}---\left(17\right);$

$\begin{array}{c}{\mathrm{\&delta;x}}_a\left(t_k\right)=\frac{1}{{n_b}^7}\{n_b\&lsqb;120{\mathrm{\&lambda;}}_{x5}{t}_k+{\mathrm{\&lambda;}}_{x4}(24-12n_b^2{t}_k^2+n_b^4{t}_k^4)+n_b^2(-2{\mathrm{\&lambda;}}_{x2}-6{\mathrm{\&lambda;}}_{x3}{t}_k+n_b^2({\mathrm{\&lambda;}}_{x0}+\\ {\mathrm{\&lambda;}}_{x1}{t}_k\left)\right)+n_b^2{t}_k^2(-20{\mathrm{\&lambda;}}_{x5}{t}_k+n_b^2({\mathrm{\&lambda;}}_{x2}+{\mathrm{\&lambda;}}_{x3}{t}_k+{\mathrm{\&lambda;}}_{x5}{t}_k^3\left)\right)\&rsqb;-n_b(24{\mathrm{\&lambda;}}_{x4}-2{\mathrm{\&lambda;}}_{x2}{n}_b^2+\\ {\mathrm{\&lambda;}}_{x0}{n}_b^4)\cos \left(n_b{t}_k\right)-(120{\mathrm{\&lambda;}}_{x5}-6{\mathrm{\&lambda;}}_{x3}{n}_b^2+{\mathrm{\&lambda;}}_{x1}{n}_b^4)\&CenterDot;\sin \left(n_b{t}_k\right)\}\end{array}---\left(18\right);$

$\begin{array}{c}{\mathrm{\&delta;y}}_a\left(t_k\right)=\frac{1}{2{n_b}^2}{\mathrm{\&lambda;}}_{y0}\left(\cosh \right(\sqrt{2}{n}_b{t}_k)-1)+\frac{1}{4{n_b}^3}{\mathrm{\&lambda;}}_{y1}(-2n_b{t}_k+\sqrt{2}\sinh (\sqrt{2}{n}_b{t}_k\left)\right)-\frac{1}{2{n_b}^4}{\mathrm{\&lambda;}}_{y2}(-\\ 1-n_b^2{t}_k^2+\cosh \left(\right(\sqrt{2}{n}_b{t}_k\left)\right)-\frac{1}{4{n_b}^5}{\mathrm{\&lambda;}}_{y3}(6n_b{t}_k+2n_b^3{t}_k^3-3\sqrt{2}\sinh (\sqrt{2}{n}_b{t}_k\left)\right)-\\ \frac{1}{2{n_b}^6}{\mathrm{\&lambda;}}_{y4}(6+6n_b^2{t}_k^2+n_b^4{t}_k^4-6\cosh (\sqrt{2}{n}_b{t}_k\left)\right)-\frac{1}{2{n_b}^7}{\mathrm{\&lambda;}}_{y5}(30n_b{t}_k+10n_b^3{t}_k^3+\\ n_b^5{t}_k^5-15\sqrt{2}\sinh \left(\sqrt{2}{n}_b{t}_k\right))\end{array}---\left(19\right);$

$\begin{array}{c}{\mathrm{\&delta;z}}_a\left(t_k\right)=\frac{1}{{n_b}^7}\{n_b\&lsqb;120{\mathrm{\&lambda;}}_{z5}{t}_k+{\mathrm{\&lambda;}}_{z4}(24-12n_b^2{t}_k^2+n_b^4{t}_k^4)+n_b^2(-2{\mathrm{\&lambda;}}_{z2}-6{\mathrm{\&lambda;}}_{z3}{t}_k+\\ n_b^2({\mathrm{\&lambda;}}_{z0}+{\mathrm{\&lambda;}}_{z1}{t}_k))+n_b^2{t}_k^2(-20{\mathrm{\&lambda;}}_{z5}{t}_k+n_b^2({\mathrm{\&lambda;}}_{z2}+{\mathrm{\&lambda;}}_{z3}{t}_k+{\mathrm{\&lambda;}}_{z5}{t}_k^3\left)\right)\&rsqb;-\\ n_b(24{\mathrm{\&lambda;}}_{z4}-2{\mathrm{\&lambda;}}_{z2}{n}_b^2+{\mathrm{\&lambda;}}_{z0}{n}_b^4)\cos \left(n_b{t}_k\right)-(120{\mathrm{\&lambda;}}_{z5}-6{\mathrm{\&lambda;}}_{z3}{n}_b^2+{\mathrm{\&lambda;}}_{z1}{n}_b^4)\&CenterDot;\sin \left(n_b{t}_k\right)\}\end{array}---\left(20\right).$

以上技术方案中优选的，所述观测向量具体是基于弹道助推段n个离散时间点对应的位置矢量求解得到的n组发射坐标系中扰动引力值，所述离散时间点按照前密后松的原则进行选取。

以上技术方案中优选的，所述第四步中视加速度偏差量求解具体为：

视加速度由气动力和发动机推力两部分组成，故而视角加速度偏差可表示为表达式(21)，详情如下：

$\mathrm{\&delta;}\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{W}}=\mathrm{\&delta;}\stackrel{\&OverBar;}{R}+\mathrm{\&delta;}\stackrel{\&OverBar;}{T}---\left(21\right);$

其中和各自的表达式详见(22)，详情如下：

$\begin{array}{c}\stackrel{\&OverBar;}{R}=\frac{1}{m}\left[\begin{array}{c}-C_x\frac{1}{2}{\mathrm{\&rho;v}}^2{S}_m\\ C_y\frac{1}{2}{\mathrm{\&rho;v}}^2{S}_m\\ C_z\frac{1}{2}{\mathrm{\&rho;v}}^2{S}_m\end{array}\right]\\ \stackrel{\&OverBar;}{T}=\frac{1}{m}\left[\begin{array}{c}-\stackrel{\&CenterDot;}{m}{u}_e+S_e(p_e-p_H)\\ 0\\ 0\end{array}\right]\end{array}---\left(22\right);$

记M_v为对速度矢量的偏导数，M_r为对位置矢量的偏导数，具体表示如表达式(23)，详情如下：

$\begin{array}{c}{M}_v=\frac{\&part;\stackrel{\&OverBar;}{R}}{\&part;{\stackrel{\&OverBar;}{V}}^T}=\frac{S_m}{m}\left[\begin{array}{ccc}-C_x{\mathrm{\&rho;v}}_x& -C_x{\mathrm{\&rho;v}}_y& -C_x{\mathrm{\&rho;v}}_z\\ C_y^{\mathrm{\&alpha;}}\mathrm{\&rho;}({\mathrm{\&alpha;v}}_x+\frac{v^2{V}_y}{2(V_x^2+V_y^2)})& C_y^{\mathrm{\&alpha;}}\mathrm{\&rho;}({\mathrm{\&alpha;v}}_y-\frac{v^2{V}_x}{2(V_x^2+V_y^2)})& C_y^{\mathrm{\&alpha;}}{\mathrm{\&rho;\&alpha;v}}_z\\ C_z^{\mathrm{\&beta;}}\mathrm{\&rho;}({\mathrm{\&beta;v}}_x-\frac{V_x{V}_z}{2\sqrt{V_x^2+V_y^2}})& C_z^{\mathrm{\&beta;}}\mathrm{\&rho;}({\mathrm{\&beta;v}}_y-\frac{V_y{V}_z}{2\sqrt{V_x^2+V_y^2}})& C_z^{\mathrm{\&beta;}}\mathrm{\&rho;}({\mathrm{\&beta;v}}_z+\frac{v^2}{2\sqrt{V_x^2+V_y^2}})\end{array}\right]\\ M_r=\frac{\&part;\stackrel{\&OverBar;}{R}}{\&part;{\stackrel{\&OverBar;}{\mathrm{\&rho;}}}^T}=\frac{S_m}{m}\left[\begin{array}{ccc}0& -C_x\frac{1}{2}{v}^2\frac{\&part;\mathrm{\&rho;}}{\&part;y}& 0\\ 0& C_y^{\mathrm{\&alpha;}}\mathrm{\&alpha;}\frac{1}{2}{v}^2\frac{\&part;\mathrm{\&rho;}}{\&part;y}& 0\\ 0& C_z^{\mathrm{\&beta;}}\mathrm{\&beta;}\frac{1}{2}{v}^2\frac{\&part;\mathrm{\&rho;}}{\&part;y}& 0\end{array}\right]\end{array}---\left(23\right);$

其中，为大气密度对高度的偏导数；

记N_r为对位置矢量的偏导数，具体详见表达式(24)：

$N_{\mathrm{\&gamma;}}=\frac{\&part;\stackrel{\&OverBar;}{T}}{\&part;{\stackrel{\&OverBar;}{\mathrm{\&rho;}}}^T}=\frac{S_e}{m}\left[\begin{array}{ccc}0& -\frac{\&part;p_H}{\&part;y}& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]---\left(24\right);$

其中，为大气压强对高度的偏导数；

弹体系到发射系、速度系到发射系以及发射惯性系到发射系的转换矩阵分别为表达式(25)、(26)以及(27)：

$G_v=\left[\begin{array}{ccc}\cos \left(\mathrm{\&theta;}\right)\cos \left(\mathrm{\&sigma;}\right)& -\sin \left(\mathrm{\&theta;}\right)& \cos \left(\mathrm{\&theta;}\right)\sin \left(\mathrm{\&sigma;}\right)\\ \sin \left(\mathrm{\&theta;}\right)\cos \left(\mathrm{\&sigma;}\right)& \cos \left(\mathrm{\&theta;}\right)& \sin \left(\mathrm{\&theta;}\right)\sin \left(\mathrm{\&sigma;}\right)\\ -\sin \left(\mathrm{\&theta;}\right)& 0& \cos \left(\mathrm{\&sigma;}\right)\end{array}\right]---\left(26\right);$

$G_A=\left[\begin{array}{ccc}1& {\mathrm{\&omega;}}_{ez}t& -{\mathrm{\&omega;}}_{ey}t\\ -{\mathrm{\&omega;}}_{ex}t& 1& {\mathrm{\&omega;}}_{ez}t\\ {\mathrm{\&omega;}}_{ey}t& -{\mathrm{\&omega;}}_{ex}t& 1\end{array}\right]---\left(27\right);$

其中，G_b中忽略了滚转角，G_v中忽略了倾侧角，G_A中各项近似至ω_et的一次项；

将表达式(23)、(24)代入表达式(21)中得到助推段弹道每一点视加速度偏差与该点的速度、位置矢量偏差之间的表达式(28)：

$\mathrm{\&delta;}\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{W}}=\&lsqb;\begin{array}{cc}{M}_v& M_r+N_r\end{array}\&rsqb;\left(\begin{array}{c}\mathrm{\&delta;}\stackrel{\&OverBar;}{V}\\ \mathrm{\&delta;}\stackrel{\&OverBar;}{\mathrm{\&rho;}}\end{array}\right)---\left(28\right).$

以上技术方案中优选的，所述第五步中解析模型迭代修正具体是：

通过不断迭代修正扰动量使得位置、速度状态矢量偏差值逐渐逼近真实值，具体迭代过程如下：

迭代初值为：取扰动量即取0；求解扰动引力并保存到数组中，迭代过程中不需重新计算；同时取

步骤一：对扰动量进行拟合；

步骤二：由表达式(15)、(16)、(17)、(18)、(19)以及(20)计算并输出从开始到发动机关机整个过程中每个时间间隔的状态偏差 ${\left[\begin{array}{cc}\mathrm{\&delta;}\stackrel{\&OverBar;}{V}\left(t\right)& \mathrm{\&delta;}\stackrel{\&OverBar;}{\mathrm{\&rho;}}\left(t\right)\end{array}\right]}^T,$ 其中输出间隔与扰动量拟合取样间隔一致；

步骤三：由表达式(28)计算步骤二中每一组输出量对应的视加速度偏差

步骤四：更新每一点扰动量的值，即

步骤五：判断 $\left\{\begin{array}{c}\left|\right|\mathrm{\&delta;}\stackrel{\&OverBar;}{V}(t_k,n)-\mathrm{\&delta;}\stackrel{\&OverBar;}{V}(t_k,n-1)\left|\right|<=eps\_v\\ \left|\right|\mathrm{\&delta;}\stackrel{\&OverBar;}{\mathrm{\&rho;}}(t_k,n)-\mathrm{\&delta;}\stackrel{\&OverBar;}{\mathrm{\&rho;}}(t_k,n-1)\left|\right|<=eps\_r\end{array}\right.$ 是否成立，即相邻两次关机点状态偏差是否小于一定值，若成立，则迭代结束；若不成立，则返回步骤一。

应用本发明的技术方案，本发明扰动引力作用下弹道助推段状态偏差解析计算方法具有“计算速度快、适应区域广、修正模型精”的特征，能满足快速机动发射及弹上实时计算的要求，具体是：

(1)本发明基于小偏差假设构建了状态偏差计算模型，模型简化合理，易于导出状态转移矩阵的解析解，为进一步推导状态偏差解析表达式打下了基础。

(2)本发明考虑了扰动引力与视加速度的耦合效应，推导了扰动引力作用下助推段状态偏差与视加速度偏差之间的映射关系，为进一步修正助推段状态偏差奠定了基础。

(3)本发明针对视加速度与状态偏差量之间转换矩阵的复杂特性，巧妙的将视加速度偏差作为扰动引力的高阶修正项，并用牛顿迭代法改善助推段状态偏差的计算精度。

(4)本发明具有计算速度快、适应区域广、修正模型精、弹上存储量小等的特征，状态偏差的计算时间控制在3s以内，解析模型计算误差与弹道求差法相比较能控制在8％以内，且计算方法适应于任意弹道。

除了上面所描述的目的、特征和优点之外，本发明还有其它的目的、特征和优点。下面将参照图，对本发明作进一步详细的说明。

附图说明

构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1为本发明算法结果与弹道求差法结果对比图；

图1a为x方向速度偏差图；

图1b为y方向速度偏差图；

图1c为z方向速度偏差图；

图1d为x方向位置偏差图；

图1e为y方向位置偏差图；

图1f为z方向位置偏差图；

图2为本发明助推段扰动量拟合精度示意图；

图2a为x方向扰动量拟合曲线图；

图2b为y方向扰动量拟合曲线图；

图2c为z方向扰动量拟合曲线图；

图2d为x方向扰动量拟合残差图；

图2e为y方向扰动量拟合残差图；

图2f为z方向扰动量拟合残差图；

图3为本发明完整解析解求解助推段速度偏差精度对比图；

图3a为x方向速度偏差图；

图3b为y方向速度偏差图；

图3c为z方向速度偏差图；

图3d为x方向相对速度误差图；

图3e为y方向相对速度误差图；

图3f为z方向相对速度误差图；

图4为本发明完整解析解求解助推段位置偏差精度对比图；

图4a为x方向位置偏差图；

图4b为y方向位置偏差图；

图4c为z方向位置偏差图；

图4d为x方向相对位置误差图；

图4e为y方向相对位置误差图；

图4f为z方向相对位置误差图；

图5为本发明考虑视加速度偏差后ΔVx的迭代效果图；

图5a为x方向速度偏差迭代效果图；

图5b为x方向速度残差修正前后对比图；

图6为本发明考虑视加速度偏差后ΔVy的迭代效果图；

图6a为y方向速度偏差迭代效果图；

图6b为y方向速度残差修正前后对比图；

图7为本发明考虑视加速度偏差后ΔVz的迭代效果图；

图7a为z方向速度偏差迭代效果图；

图7b为z方向速度残差修正前后对比图。

具体实施方式

以下结合附图对本发明的实施例进行详细说明，但是本发明可以根据权利要求限定和覆盖的多种不同方式实施。

实施例1：

以某型号弹道导弹为仿真对象，仿真初始条件设置如表1所示：

表1仿真发射点的参数值

发射点经度(°)	发射点纬度(°)	发射点高度(m)	发射方位角(°)	关机点时间(s)
					86.5E	31.5N	5584	30.0	176.5

扰动计算方法为：球谐函数法。

引力模型为：截断到360阶的EGM2008模型。

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

具体包含以下步骤(各表达式中的参数所代表的含义如表3所示)：

第一步，助推段标准弹道设计，具体是：

按照发射点、目标点、某一型号导弹总体参数、大气模型以及地球引力模型等任务条件，设计一条满足条件的弹道，并按照一定时间间隔(一般取2s)保存助推段发射系中时间、速度、位置、质量、程序角等状态值，这里时间间隔取2s，需要保存89组数据。

第二步，状态偏差模型建立，具体是：

标准弹道计算模型可以表示为表达式(1)，详情如下：

$\left\{\begin{array}{c}{\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{V}}}^{*}={\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{W}}}^{*}+g^{*}\left({\stackrel{\&OverBar;}{\mathrm{\&rho;}}}^{*}\right)\\ {\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{\mathrm{\&rho;}}}}^{*}={\stackrel{\&OverBar;}{V}}^{*}\end{array}\right.---\left(1\right);$

其中，为标准弹道视加速度，为标准弹道引力加速度；

若考虑扰动引力对导弹运动的影响，则实际飞行弹道与标准弹道的等时变分可表示为表达式(2)，详情如下：

$\left\{\begin{array}{c}\mathrm{\&delta;}\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{V}}=\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{W}}-{\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{W}}}^{*}+\stackrel{\&OverBar;}{g}\left(\stackrel{\&OverBar;}{\mathrm{\&rho;}}\right)-{\stackrel{\&OverBar;}{g}}^{*}\left({\stackrel{\&OverBar;}{\mathrm{\&rho;}}}^{*}\right)\\ \mathrm{\&delta;}\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{\mathrm{\&rho;}}}=\mathrm{\&delta;}\stackrel{\&OverBar;}{V}\end{array}\right.---\left(2\right);$

其中，为实际弹道视加速度，为实际弹道引力加速度。

将表达式(2)改写为矩阵形式，同时经过小偏差处理并略去高阶小量(具体为：采用常规非线性系统线性化的处理方法，即泰勒展开后取一次项，忽略高阶项)，可得表达式(3)，详情如下：

$\left(\begin{array}{c}\mathrm{\&delta;}\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{V}}\\ \mathrm{\&delta;}\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{\mathrm{\&rho;}}}\end{array}\right)=T\left(\begin{array}{c}\mathrm{\&delta;}\stackrel{\&OverBar;}{V}\\ \mathrm{\&delta;}\stackrel{\&OverBar;}{\mathrm{\&rho;}}\end{array}\right)+\left(\begin{array}{c}\mathrm{\&delta;}\stackrel{\&OverBar;}{F}\\ 0\end{array}\right)---\left(3\right);$

其中，为扰动量，且有为扰动引力，为视加速度偏差(视加速度偏差由两部分组成，一部分为惯性工具误差导致的视加速度偏差；另一部分为扰动引力与视加速度耦合作用引起的视加速度偏差，本实施例只分析扰动引力对弹道助推段状态偏差的影响，因而忽略惯性工具误差导致的视加速度偏差)；式中T的表达式详见(4)，详情如下：

表达式(4)中A₁、A₂、A₃、A₄、A₅以及A₆六者各自的表达式详见(5)：

$\left\{\begin{array}{c}{A}_1=-n_b^2(1-3\frac{x^2}{r^2})\\ A_2=3n_b^2\frac{x}{r}\frac{R_0+y}{r}\\ A_3=3n_b^2\frac{x}{r}\&CenterDot;\frac{z}{r}\\ A_4=-n_b^2(1-3\frac{{(R_0+y)}^2}{r^2})\\ A_5=3n_b^2\frac{z}{r}\&CenterDot;\frac{R_0+y}{r}\\ A_6=-n_b^2(1-3\frac{z^2}{r^2})\end{array}\right.---\left(5\right);$

经过合理简化，采用伴随方程求得表达式(3)的状态转移矩阵解析解为表达式(7)，详情如下：

$\begin{array}{c}\mathrm{\&Phi;}(t_k,t)=\mathrm{\&Phi;}\left(\mathrm{\&tau;}\right)\\ =\left[\begin{array}{cccccc}{\mathrm{cosn}}_b\mathrm{\&tau;}& {\mathrm{\&Phi;}}_{12}\left(\mathrm{\&tau;}\right)& 0& -n_b\sin n_b\mathrm{\&tau;}& 0& 0\\ {\mathrm{\&Phi;}}_{21}\left(\mathrm{\&tau;}\right)& \cosh \sqrt{2}{n}_b\mathrm{\&tau;}& 0& 0& \sqrt{2}{n}_b\sinh \sqrt{2}{n}_b\mathrm{\&tau;}& 0\\ 0& 0& \cos n_b\mathrm{\&tau;}& 0& 0& -n_b\sin n_b\mathrm{\&tau;}\\ \frac{\sin n_b\mathrm{\&tau;}}{n_b}& 0& 0& {\mathrm{cosn}}_b\mathrm{\&tau;}& 0& 0\\ 0& \frac{\sinh \sqrt{2}{n}_b\mathrm{\&tau;}}{\sqrt{2}{n}_b}& 0& 0& \cosh \sqrt{2}{n}_b\mathrm{\&tau;}& 0\\ 0& 0& \frac{\sin n_b\mathrm{\&tau;}}{n_b}& 0& 0& \cos n_b\mathrm{\&tau;}\end{array}\right]\end{array}---\left(7\right);$

其中Φ₁₂(τ)的表达式详见(8)：

${\mathrm{\&Phi;}}_{12}\left(\mathrm{\&tau;}\right)={\mathrm{\&Phi;}}_{21}\left(\mathrm{\&tau;}\right)=\frac{{\left(n_b\mathrm{\&tau;}\right)}^2}{3!}(1+\frac{2}{3}\frac{\mathrm{\&tau;}}{t_k}+\frac{1}{6}\frac{{\mathrm{\&tau;}}^2}{t_k^2})---\left(8\right).$

所采用的伴随方程的具体解法如下：

定义伴随方程为： $\frac{d}{dt}G(t,\mathrm{\&tau;})=-T^T\left(t\right)G(t,\mathrm{\&tau;});$

由于则G^T(t,τ)＝Φ^-1(t,τ)＝Φ(τ，t)；

积分方程 $\left\{\begin{array}{c}\frac{dG(t,t_k)}{dt}=-T^T\left(t\right)G(t,t_k)\\ G(t_k,t_k)=I_6\end{array}\right.$ 从t_k到t₀积分，即可解得G(t,t_k)，则Φ(t_k,t)＝G^T(t,t_k)。

上述合理简化具体过程如下：

由表达式(5)可见，A₁、A₂、A₃、A₄、A₅以及A₆六者均为小量，其中A₁、A₄以及A₆较大，其数量级与相当，而A₂比小一个数量级，A₃和A₅则比小两三个数量级，即有：

$\left\{\begin{array}{c}\frac{R_0+y}{r}\&ap;1\\ \frac{z^2}{r^3}<<1\\ \frac{x^2}{r^2}<<1\end{array}\right.$

因此可将表达式(5)改写为表达式(6)的形式，详情如下：

$\left\{\begin{array}{c}{A}_1=A_6=-n_b^2\\ A_4=2n_b^2\\ A_2=\frac{3x}{r}{n}_b^2\\ A_3=A_5=0\end{array}\right.---\left(6\right);$

同时，考虑到n_b在整个助推段过程中的变化范围小于1％，因此积分时设n_b为常数。

根据表达式(7)即可推得助推段扰动引力引起的状态偏差半解析表达式(9)，如下：

$\begin{array}{c}{\mathrm{\&delta;v}}_{ax}\left(t_k\right)={\&Integral;}_0^{t_k}{\mathrm{\&Phi;}}_{11}(t_k-\mathrm{\&tau;}){\mathrm{\&delta;F}}_x\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;}+{\&Integral;}_0^{t_k}{\mathrm{\&Phi;}}_{12}(t_k-\mathrm{\&tau;}){\mathrm{\&delta;F}}_y\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;}\\ {\mathrm{\&delta;v}}_{ay}\left(t_k\right)={\&Integral;}_0^{t_k}{\mathrm{\&Phi;}}_{22}(t_k-\mathrm{\&tau;}){\mathrm{\&delta;F}}_y\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;}+{\&Integral;}_0^{t_k}{\mathrm{\&Phi;}}_{21}(t_k-\mathrm{\&tau;}){\mathrm{\&delta;F}}_x\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;}\\ {\mathrm{\&delta;v}}_{az}\left(t_k\right)={\&Integral;}_0^{t_k}{\mathrm{\&Phi;}}_{33}(t_k-\mathrm{\&tau;}){\mathrm{\&delta;F}}_z\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;}\\ {\mathrm{\&delta;x}}_a\left(t_k\right)={\&Integral;}_0^{t_k}{\mathrm{\&Phi;}}_{41}(t_k-\mathrm{\&tau;}){\mathrm{\&delta;F}}_x\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;}\\ {\mathrm{\&delta;y}}_a\left(t_k\right)={\&Integral;}_0^{t_k}{\mathrm{\&Phi;}}_{52}(t_k-\mathrm{\&tau;}){\mathrm{\&delta;F}}_y\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;}\\ {\mathrm{\&delta;z}}_a\left(t_k\right)={\&Integral;}_0^{t_k}{\mathrm{\&Phi;}}_{63}(t_k-\mathrm{\&tau;}){\mathrm{\&delta;F}}_z\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;}\end{array}---\left(9\right).$

为了验证表达式(9)的计算精度，将其计算(简称“半解析法”，下文同)结果与弹道求差法(简称“求差法”)的结果进行对比详见图1(图1中，图1a为x方向速度偏差图，图1b为y方向速度偏差图,图1c为z方向速度偏差图,图1d为x方向位置偏差图，图1e为y方向位置偏差图,图1f为z方向位置偏差图)，其中，弹道求差法具体是：采用四阶龙格-库塔法积分求解两次弹道，其中两条弹道计算中引力模型一条采用正常引力模型(只考虑到J2项)，另一条考虑360阶扰动引力，同时对弹道每一点积分状态值等时求差，即为扰动引力作用下弹道导弹助推段状态偏差值。从图1中可以看出，本发明提出的半解析计算方法基本能反映扰动引力对助推段状态偏差的影响特性，但由于忽略了视加速度与扰动引力的耦合影响，半解析法所求结果相对求差法来说误差较大，部分时间点对应的速度矢量相对误差超过20％，因而不能满足实际助推段状态偏差精确计算的要求。

第三步、助推段扰动引力拟合，具体为：

采用最小二乘法将助推段扰动引力拟合为关于时间的多项式函数，且发射惯性系中扰动引力三分量分别进行拟合，具体是：设观测方程为表达式(10)：

Z_n×1＝H_n×(m+1)λ_(m+1)×1+ε(10)；

其中，Z_n×1为观测向量(观测向量即为基于弹道助推段n个离散时间点对应的位置矢量求解得到的n组发射坐标系中扰动引力值，离散时间点的选取按照前密后松的原则，即越接近地面扰动引力波动越大，相应的离散点应该越密)，n为拟合点的个数；H_n×(m+1)为系数矩阵，m代表了多项式拟合次数；λ_(m+1)×1为待定参数向量；ε为随机误差向量。

拟合多项式采用经验公式(11)则可得系数矩阵以及扰动引力拟合系数分别为表述式(12)和(13)，详情如下：

$\mathrm{\&delta;}F\left(t\right)=\underset{i=0}{\overset{n}{\&Sigma;}}{\mathrm{\&lambda;}}_i{\left(\frac{t}{10}\right)}^i---\left(11\right);$

λ＝(H^TH)^-1H^TZ(13)。

理论上拟合多项式的阶次越高拟合精度越高，实际上当次数高于一定值时便会产生Runge现象。对于本实施例中扰动引力的拟合，拟合次数取5阶即可满足较高精度(本实施例的拟合精度详见图2所示(图2中，图2a为x方向扰动量拟合曲线图，图2b为y方向扰动量拟合曲线图,图2c为z方向扰动量拟合曲线图,图2d为x方向扰动量拟合残差图，图2e为y方向扰动量拟合残差图,图2f为z方向扰动量拟合残差图)，从图2可知，拟合精度在0.01mgal量级，扰动引力的高精度拟合为解析表达式的高精度奠定了基础)，若要获得更高精度，最好采用分段拟合而不是增加多项式拟合阶次。

当拟合多项式取5阶时，可得扰动量拟合表达式(14)，如下：

$\left\{\begin{array}{c}{\mathrm{\&delta;F}}_x\left(t\right)={\mathrm{\&lambda;}}_{x0}+{\mathrm{\&lambda;}}_{x1}(t/10)+{\mathrm{\&lambda;}}_{x2}{(t/10)}^2+{\mathrm{\&lambda;}}_{x3}{(t/10)}^3+{\mathrm{\&lambda;}}_{x4}{(t/10)}^4+{\mathrm{\&lambda;}}_{x5}{(t/10)}^5\\ {\mathrm{\&delta;F}}_y\left(t\right)={\mathrm{\&lambda;}}_{y0}+{\mathrm{\&lambda;}}_{y1}(t/10)+{\mathrm{\&lambda;}}_{y2}{(t/10)}^2+{\mathrm{\&lambda;}}_{y3}{(t/10)}^3+{\mathrm{\&lambda;}}_{y4}{(t/10)}^4+{\mathrm{\&lambda;}}_{y5}{(t/10)}^5\\ {\mathrm{\&delta;F}}_z\left(t\right)={\mathrm{\&lambda;}}_{z0}+{\mathrm{\&lambda;}}_{z1}(t/10)+{\mathrm{\&lambda;}}_{z2}{(t/10)}^2+{\mathrm{\&lambda;}}_{z3}{(t/10)}^3+{\mathrm{\&lambda;}}_{z4}{(t/10)}^4+{\mathrm{\&lambda;}}_{z5}{(t/10)}^5\end{array}\right.---\left(14\right).$

将表达式(14)带入助推段扰动引力引起的状态偏差半解析表达式(9)导出助推段扰动引力引起的状态偏差完整解析表达式(15)、(16)、(17)、(18)、(19)以及(20)，详情如下：

$\begin{array}{c}{\mathrm{\&delta;v}}_{ax}\left(t_k\right)=\frac{1}{n_b}{\mathrm{\&lambda;}}_{x0}\sin \left(n_b{t}_k\right)+\frac{1}{{n_b}^2}{\mathrm{\&lambda;}}_{x1}(1-\cos (n_b{t}_k)+\frac{2}{{n_b}^3}{\mathrm{\&lambda;}}_{x2}(n_b{t}_k-\sin \left(n_b{t}_k\right))+\\ \frac{3}{{n_b}^4}{\mathrm{\&lambda;}}_{x3}(-2+n_b^2{t}_k^2+2\cos (n_b{t}_k\left)\right)+\frac{4}{{n_b}^5}{\mathrm{\&lambda;}}_{x4}(-6n_b{t}_k+n_b^3{t}_k^3+6\sin (n_b{t}_k\left)\right)+\\ \frac{5}{{n_b}^6}{\mathrm{\&lambda;}}_{x5}(24-12n_b^2{t}_k^2+n_b^4{t}_k^4-24\cos (n_b{t}_k\left)\right)+\frac{1}{90720}\&CenterDot;\{n_b^2{t}_k^3\&CenterDot;\&lsqb;8064{\mathrm{\&lambda;}}_{y0}+\\ t_k\&CenterDot;(1848{\mathrm{\&lambda;}}_{y1}+(696{\mathrm{\&lambda;}}_{y2}+333{\mathrm{\&lambda;}}_{y3}{t}_k+184333{\mathrm{\&lambda;}}_{y4}{t}_k^2+112{\mathrm{\&lambda;}}_{y5}{t}_k^3\left)\right)\&rsqb;\}\end{array}---\left(15\right);$

$\begin{array}{c}{\mathrm{\&delta;v}}_{ay}\left(t_k\right)=\frac{1}{\sqrt{2}{n}_b}{\mathrm{\&lambda;}}_{y0}\sinh \left(\sqrt{2}{n}_b{t}_k\right)+\frac{1}{2{n_b}^2}{\mathrm{\&lambda;}}_{y1}\left(\cosh \right(\sqrt{2}{n}_b{t}_k)-1)+\frac{2}{2{n_b}^3}{\mathrm{\&lambda;}}_{y2}(-2n_b{t}_k+\\ \sqrt{2}\sinh \left(\sqrt{2}{n}_b{t}_k\right))+\frac{3}{2{n_b}^4}{\mathrm{\&lambda;}}_{y3}(-1-n_b^2{t}_k^2+\cosh (\sqrt{2}{n}_b{t}_k\left)\right)+\frac{1}{{n_b}^5}{\mathrm{\&lambda;}}_{y4}(-2n_b{t}_k(3+\\ n_b^2{t}_k^2)+3\sqrt{2}\sinh \left(\right(\sqrt{2}{n}_b{t}_k\left)\right)-\frac{5}{{n_b}^6}{\mathrm{\&lambda;}}_{y5}(6+12n_b^2{t}_k^2+n_b^4{t}_k^4-6\cosh (\sqrt{2}{n}_b{t}_k\left)\right)+\\ \frac{1}{90720}\&CenterDot;\{n_b^2{t}_k^3\&CenterDot;\&lsqb;8064{\mathrm{\&lambda;}}_{x0}+t_k\&CenterDot;(1848{\mathrm{\&lambda;}}_{x1}+(696{\mathrm{\&lambda;}}_{x2}+333{\mathrm{\&lambda;}}_{x3}{t}_k+184333{\mathrm{\&lambda;}}_{x4}{t}_k^2+\\ 112{\mathrm{\&lambda;}}_{x5}{t}_k^3\left)\right)\&rsqb;\}\end{array}---\left(16\right);$

$\begin{array}{c}{\mathrm{\&delta;v}}_{az}\left(t_k\right)=\frac{1}{n_b}{\mathrm{\&lambda;}}_{z0}\sin \left(n_b{t}_k\right)+\frac{1}{{n_b}^2}{\mathrm{\&lambda;}}_{z1}(1-\cos (n_b{t}_k)+\frac{2}{{n_b}^3}{\mathrm{\&lambda;}}_{z2}(n_b{t}_k-\sin \left(n_b{t}_k\right))+\\ \frac{3}{{n_b}^4}{\mathrm{\&lambda;}}_{z3}(-2+n_b^2{t}_k^2+2\cos (n_b{t}_k\left)\right)+\frac{4}{{n_b}^5}{\mathrm{\&lambda;}}_{z4}(-6n_b{t}_k+n_b^3{t}_k^3+6\sin (n_b{t}_k\left)\right)+\\ \frac{5}{{n_b}^6}{\mathrm{\&lambda;}}_{z5}(24-12n_b^2{t}_k^2+n_b^4{t}_k^4-24\cos (n_b{t}_k\left)\right)\end{array}---\left(17\right);$

$\begin{array}{c}{\mathrm{\&delta;x}}_a\left(t_k\right)=\frac{1}{{n_b}^7}\{n_b\&lsqb;120{\mathrm{\&lambda;}}_{x5}{t}_k+{\mathrm{\&lambda;}}_{x4}(24-12n_b^2{t}_k^2+n_b^4{t}_k^4)+n_b^2(-2{\mathrm{\&lambda;}}_{x2}-6{\mathrm{\&lambda;}}_{x3}{t}_k+n_b^2({\mathrm{\&lambda;}}_{x0}+\\ {\mathrm{\&lambda;}}_{x1}{t}_k\left)\right)+n_b^2{t}_k^2(-20{\mathrm{\&lambda;}}_{x5}{t}_k+n_b^2({\mathrm{\&lambda;}}_{x2}+{\mathrm{\&lambda;}}_{x3}{t}_k+{\mathrm{\&lambda;}}_{x5}{t}_k^3\left)\right)\&rsqb;-n_b(24{\mathrm{\&lambda;}}_{x4}-2{\mathrm{\&lambda;}}_{x2}{n}_b^2+\\ {\mathrm{\&lambda;}}_{x0}{n}_b^4)\cos \left(n_b{t}_k\right)-(120{\mathrm{\&lambda;}}_{x5}-6{\mathrm{\&lambda;}}_{x3}{n}_b^2+{\mathrm{\&lambda;}}_{x1}{n}_b^4)\&CenterDot;\sin \left(n_b{t}_k\right)\}\end{array}---\left(18\right);$

$\begin{array}{c}{\mathrm{\&delta;y}}_a\left(t_k\right)=\frac{1}{2{n_b}^2}{\mathrm{\&lambda;}}_{y0}\left(\cosh \right(\sqrt{2}{n}_b{t}_k)-1)+\frac{1}{4{n_b}^3}{\mathrm{\&lambda;}}_{y1}(-2n_b{t}_k+\sqrt{2}\sinh (\sqrt{2}{n}_b{t}_k\left)\right)-\frac{1}{2{n_b}^4}{\mathrm{\&lambda;}}_{y2}(-\\ 1-n_b^2{t}_k^2+\cosh \left(\right(\sqrt{2}{n}_b{t}_k\left)\right)-\frac{1}{4{n_b}^5}{\mathrm{\&lambda;}}_{y3}(6n_b{t}_k+2n_b^3{t}_k^3-3\sqrt{2}\sinh (\sqrt{2}{n}_b{t}_k\left)\right)-\\ \frac{1}{2{n_b}^6}{\mathrm{\&lambda;}}_{y4}(6+6n_b^2{t}_k^2+n_b^4{t}_k^4-6\cosh (\sqrt{2}{n}_b{t}_k\left)\right)-\frac{1}{2{n_b}^7}{\mathrm{\&lambda;}}_{y5}(30n_b{t}_k+10n_b^3{t}_k^3+\\ n_b^5{t}_k^5-15\sqrt{2}\sinh \left(\sqrt{2}{n}_b{t}_k\right))\end{array}---\left(19\right);$

$\begin{array}{c}{\mathrm{\&delta;z}}_a\left(t_k\right)=\frac{1}{{n_b}^7}\{n_b\&lsqb;120{\mathrm{\&lambda;}}_{z5}{t}_k+{\mathrm{\&lambda;}}_{z4}(24-12n_b^2{t}_k^2+n_b^4{t}_k^4)+n_b^2(-2{\mathrm{\&lambda;}}_{z2}-6{\mathrm{\&lambda;}}_{z3}{t}_k+\\ n_b^2({\mathrm{\&lambda;}}_{z0}+{\mathrm{\&lambda;}}_{z1}{t}_k))+n_b^2{t}_k^2(-20{\mathrm{\&lambda;}}_{z5}{t}_k+n_b^2({\mathrm{\&lambda;}}_{z2}+{\mathrm{\&lambda;}}_{z3}{t}_k+{\mathrm{\&lambda;}}_{z5}{t}_k^3\left)\right)\&rsqb;-\\ n_b(24{\mathrm{\&lambda;}}_{z4}-2{\mathrm{\&lambda;}}_{z2}{n}_b^2+{\mathrm{\&lambda;}}_{z0}{n}_b^4)\cos \left(n_b{t}_k\right)-(120{\mathrm{\&lambda;}}_{z5}-6{\mathrm{\&lambda;}}_{z3}{n}_b^2+{\mathrm{\&lambda;}}_{z1}{n}_b^4)\&CenterDot;\sin \left(n_b{t}_k\right)\}\end{array}---\left(20\right).$

完整解析解的求解(简称“解析法”)精度与半解析解的求解(简称“半解析法”)精度对比如图3(图3中，图3a为x方向速度偏差图，图3b为y方向速度偏差图,图3c为z方向速度偏差图,图3d为x方向相对速度误差图，图3e为y方向相对速度误差图,图3f为z方向相对速度误差图)和图4(图4中，图4a为x方向位置偏差图，图4b为y方向位置偏差图,图4c为z方向位置偏差图,图4d为x方向相对位置误差图，图4e为y方向相对位置误差图,图4f为z方向相对位置误差图)所示，由图3和图4可知，本发明所提助推段状态偏差完整解析解相对半解析解而言精度很高，相对误差控制在2％以内，满足弹道导弹助推段状态偏差快速、精确求解的要求。

弹道求差法、状态偏差半解析方法以及完整解析方法在求解助推段状态偏差时所用时间对比如表2所示。

表2不同计算方法所用时间对比

从表2可知，半解析计算方法比弹道求差法快了4倍多，完整解析计算方法比弹道求差法快了50倍左右，因此本发明提出的方法在计算速度上具有极大的优势。

第四步、视加速度偏差量求解，具体为：

视加速度由气动力和发动机推力两部分组成，故而视角加速度偏差可表示为表达式(21)，详情如下：

$\mathrm{\&delta;}\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{W}}=\mathrm{\&delta;}\stackrel{\&OverBar;}{R}+\mathrm{\&delta;}\stackrel{\&OverBar;}{T}---\left(21\right);$

其中和各自的表达式详见(22)，详情如下：

$\begin{array}{c}\stackrel{\&OverBar;}{R}=\frac{1}{m}\left[\begin{array}{c}-C_x\frac{1}{2}{\mathrm{\&rho;v}}^2{S}_m\\ C_y\frac{1}{2}{\mathrm{\&rho;v}}^2{S}_m\\ C_z\frac{1}{2}{\mathrm{\&rho;v}}^2{S}_m\end{array}\right]\\ \stackrel{\&OverBar;}{T}=\frac{1}{m}\left[\begin{array}{c}-\stackrel{\&CenterDot;}{m}{u}_e+S_e(p_e-p_H)\\ 0\\ 0\end{array}\right]\end{array}---\left(22\right);$

记M_v为对速度矢量的偏导数，M_r为对位置矢量的偏导数，具体表示如表达式(23)，详情如下：

$\begin{array}{c}{M}_v=\frac{\&part;\stackrel{\&OverBar;}{R}}{\&part;{\stackrel{\&OverBar;}{V}}^T}=\frac{S_m}{m}\left[\begin{array}{ccc}-C_x{\mathrm{\&rho;v}}_x& -C_x{\mathrm{\&rho;v}}_y& -C_x{\mathrm{\&rho;v}}_z\\ C_y^{\mathrm{\&alpha;}}\mathrm{\&rho;}({\mathrm{\&alpha;v}}_x+\frac{v^2{V}_y}{2(V_x^2+V_y^2)})& C_y^{\mathrm{\&alpha;}}\mathrm{\&rho;}({\mathrm{\&alpha;v}}_y-\frac{v^2{V}_x}{2(V_x^2+V_y^2)})& C_y^{\mathrm{\&alpha;}}{\mathrm{\&rho;\&alpha;v}}_z\\ C_z^{\mathrm{\&beta;}}\mathrm{\&rho;}({\mathrm{\&beta;v}}_x-\frac{V_x{V}_z}{2\sqrt{V_x^2+V_y^2}})& C_z^{\mathrm{\&beta;}}\mathrm{\&rho;}({\mathrm{\&beta;v}}_y-\frac{V_y{V}_z}{2\sqrt{V_x^2+V_y^2}})& C_z^{\mathrm{\&beta;}}\mathrm{\&rho;}({\mathrm{\&beta;v}}_z+\frac{v^2}{2\sqrt{V_x^2+V_y^2}})\end{array}\right]\\ M_r=\frac{\&part;\stackrel{\&OverBar;}{R}}{\&part;{\stackrel{\&OverBar;}{\mathrm{\&rho;}}}^T}=\frac{S_m}{m}\left[\begin{array}{ccc}0& -C_x\frac{1}{2}{v}^2\frac{\&part;\mathrm{\&rho;}}{\&part;y}& 0\\ 0& C_y^{\mathrm{\&alpha;}}\mathrm{\&alpha;}\frac{1}{2}{v}^2\frac{\&part;\mathrm{\&rho;}}{\&part;y}& 0\\ 0& C_z^{\mathrm{\&beta;}}\mathrm{\&beta;}\frac{1}{2}{v}^2\frac{\&part;\mathrm{\&rho;}}{\&part;y}& 0\end{array}\right]\end{array}---\left(23\right);$

其中，为大气密度对高度的偏导数。

记N_r为对位置矢量的偏导数，具体详见表达式(24)：

$N_{\mathrm{\&gamma;}}=\frac{\&part;\stackrel{\&OverBar;}{T}}{\&part;{\stackrel{\&OverBar;}{\mathrm{\&rho;}}}^T}=\frac{S_e}{m}\left[\begin{array}{ccc}0& -\frac{\&part;p_H}{\&part;y}& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]---\left(24\right);$

其中，为大气压强对高度的偏导数。

显然，表达式(23)和(24)分别为速度系和弹体系里面的表达式，计算时需将其转换到发射坐标系中，再由发射系转换到发射惯性系中。弹体系到发射系、速度系到发射系以及发射惯性系到发射系的转换矩阵分别为(25)、(26)以及(27)：

$G_v=\left[\begin{array}{ccc}\cos \left(\mathrm{\&theta;}\right)\cos \left(\mathrm{\&sigma;}\right)& -\sin \left(\mathrm{\&theta;}\right)& \cos \left(\mathrm{\&theta;}\right)\sin \left(\mathrm{\&sigma;}\right)\\ \sin \left(\mathrm{\&theta;}\right)\cos \left(\mathrm{\&sigma;}\right)& \cos \left(\mathrm{\&theta;}\right)& \sin \left(\mathrm{\&theta;}\right)\sin \left(\mathrm{\&sigma;}\right)\\ -\sin \left(\mathrm{\&theta;}\right)& 0& \cos \left(\mathrm{\&sigma;}\right)\end{array}\right]---\left(26\right);$

$G_A=\left[\begin{array}{ccc}1& {\mathrm{\&omega;}}_{ez}t& -{\mathrm{\&omega;}}_{ey}t\\ -{\mathrm{\&omega;}}_{ex}t& 1& {\mathrm{\&omega;}}_{ez}t\\ {\mathrm{\&omega;}}_{ey}t& -{\mathrm{\&omega;}}_{ex}t& 1\end{array}\right]---\left(27\right);$

其中，G_b中忽略了滚转角，G_v中忽略了倾侧角，G_A中各项近似至ω_et的一次项。

将表达式(23)、(24)代入表达式(21)可得到助推段弹道每一点视加速度偏差与该点的速度、位置矢量偏差之间的关系表达式(28)：

$\mathrm{\&delta;}\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{W}}=\&lsqb;\begin{array}{cc}{M}_v& M_r+N_r\end{array}\&rsqb;\left(\begin{array}{c}\mathrm{\&delta;}\stackrel{\&OverBar;}{V}\\ \mathrm{\&delta;}\stackrel{\&OverBar;}{\mathrm{\&rho;}}\end{array}\right)---\left(28\right).$

需要注意的是，由于需要在发射坐标系中计算，在运用表达式(28)求解当前时刻视加速度偏差时，需将和先由发射惯性系转换到发射系中，待计算完成后还需要将再次转换到发射惯性系中。

第五步、解析模型迭代修正，具体是：

通过不断迭代修正扰动量使得位置、速度状态矢量偏差值逐渐逼近真实值，迭代效果详见图5(图5a为x方向速度偏差迭代效果图，图5b为x方向速度残差修正前后对比图)、图6(图6a为y方向速度偏差迭代效果图，图6b为y方向速度残差修正前后对比图)以及图7(图7a为z方向速度偏差迭代效果图，图7b为z方向速度残差修正前后对比图)。具体迭代过程如下：

迭代初值为：取扰动量即取0；求解扰动引力并保存到数组中，迭代过程中不需重新计算；同时取

步骤一：对扰动量进行拟合；

步骤二：由表达式(15)、(16)、(17)、(18)、(19)以及(20)计算并输出从开始到发动机关机整个过程中每个时间间隔的状态偏差 ${\left[\begin{array}{cc}\mathrm{\&delta;}\stackrel{\&OverBar;}{V}\left(t\right)& \mathrm{\&delta;}\stackrel{\&OverBar;}{\mathrm{\&rho;}}\left(t\right)\end{array}\right]}^T,$ 其中输出间隔与扰动量拟合取样间隔一致；

步骤三：由表达式(28)计算步骤二中每一组输出量对应的视加速度偏差

步骤四：更新每一点扰动量的值，即

步骤五：判断 $\left\{\begin{array}{c}\left|\right|\mathrm{\&delta;}\stackrel{\&OverBar;}{V}(t_k,n)-\mathrm{\&delta;}\stackrel{\&OverBar;}{V}(t_k,n-1)\left|\right|<=eps\_v\\ \left|\right|\mathrm{\&delta;}\stackrel{\&OverBar;}{\mathrm{\&rho;}}(t_k,n)-\mathrm{\&delta;}\stackrel{\&OverBar;}{\mathrm{\&rho;}}(t_k,n-1)\left|\right|<=eps\_r\end{array}\right.$ 是否成立，即相邻两次关机点状态偏差是否小于一定值，若成立，则迭代结束；若不成立，则返回步骤一。

从图5、图6以及图7可知，本发明所提修正方法效果明显，修正后的残差相对修正前缩小了近3倍，极大的提高了助推段状态偏差解析计算的精度。

综上所述，本发明采用的扰动引力作用下弹道助推段状态偏差解析计算方法具有以下效果：

(1)本发明基于小偏差假设构建了状态偏差计算模型，模型简化合理，易于导出状态转移矩阵的解析解，为进一步推导状态偏差解析表达式打下了基础。

(2)本发明考虑了扰动引力与视加速度的耦合效应，推导了扰动引力作用下助推段状态偏差与视加速度偏差之间的映射关系，为进一步修正助推段状态偏差奠定了基础；

(3)本发明针对视加速度与状态偏差量之间转换矩阵的复杂特性，巧妙的将视加速度偏差作为扰动引力的高阶修正项，并用牛顿迭代法改善助推段状态偏差的计算精度。

(4)本发明具有计算速度快、适应区域广、修正模型精、弹上存储量小等的特征，状态偏差的计算时间控制在3s以内，解析模型计算误差相对求差法而言控制在8％以内，且计算方法适应与任意弹道。

本实施例涉及的参数及其所具有的含义详见表3。

表3本实施涉及的参数及其含义

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

Claims (8)

1.一种扰动引力作用下弹道助推段状态偏差解析计算方法，其特征在于：包括以下步骤：第一步，助推段标准弹道设计；第二步，状态偏差模型建立；第三步，助推段扰动量拟合；第四步，视加速度偏差量求解；第五步，解析模型迭代修正。

2.根据权利要求1所述的扰动引力作用下弹道助推段状态偏差解析计算方法，其特征在于：包括以下步骤：

第一步，助推段标准弹道设计，具体为：按照发射点、目标点、某一型号导弹的总体参数、大气模型以及地球引力模型的任务条件，设计一条满足条件的弹道，并按照一定时间间隔保存助推段发射系中时间、速度、位置、质量、程序角五项状态值；

第二步，状态偏差模型建立，具体为：根据标准弹道计算模型经过一定处理获得助推段扰动引力引起的状态偏差半解析表达式(9)，详情如下：

${\mathrm{\&delta;v}}_{ax}\left(t_k\right)={\&Integral;}_0^{t_k}{\mathrm{\&Phi;}}_{11}(t_k-\mathrm{\&tau;}){\mathrm{\&delta;F}}_x\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;}+{\&Integral;}_0^{t_k}{\mathrm{\&Phi;}}_{12}(t_k-\mathrm{\&tau;}){\mathrm{\&delta;F}}_y\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;}$

${\mathrm{\&delta;v}}_{ay}\left(t_k\right)={\&Integral;}_0^{t_k}{\mathrm{\&Phi;}}_{22}(t_k-\mathrm{\&tau;}){\mathrm{\&delta;F}}_y\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;}+{\&Integral;}_0^{t_k}{\mathrm{\&Phi;}}_{21}(t_k-\mathrm{\&tau;}){\mathrm{\&delta;F}}_x\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;}$

${\mathrm{\&delta;v}}_{az}\left(t_k\right)={\&Integral;}_0^{t_k}{\mathrm{\&Phi;}}_{33}(t_k-\mathrm{\&tau;}){\mathrm{\&delta;F}}_z\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;}$

${\mathrm{\&delta;x}}_a\left(t_k\right)={\&Integral;}_0^{t_k}{\mathrm{\&Phi;}}_{41}(t_k-\mathrm{\&tau;}){\mathrm{\&delta;F}}_x\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;}---\left(9\right);$

${\mathrm{\&delta;y}}_a\left(t_k\right)={\&Integral;}_0^{t_k}{\mathrm{\&Phi;}}_{52}(t_k-\mathrm{\&tau;}){\mathrm{\&delta;F}}_y\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;}$

${\mathrm{\&delta;z}}_a\left(t_k\right)={\&Integral;}_0^{t_k}{\mathrm{\&Phi;}}_{63}(t_k-\mathrm{\&tau;}){\mathrm{\&delta;F}}_z\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;}$

第三步，助推段扰动量拟合，获得助推段扰动引力引起的状态偏差完整解析表达式，具体为：采用最小二乘法将助推段扰动量拟合为关于时间的多项式函数，且发射惯性系中扰动量三分量分别进行拟合，得到拟合多项式取5阶时的扰动量拟合表达式(14)，将扰动量拟合表达式(14)带入所述第二步中的助推段扰动引力引起的状态偏差半解析表达式(9)中获得助推段扰动引力引起的状态偏差完整解析表达式(15)、(16)、(17)、(18)、(19)以及(20)，详情如下：

$\{\begin{array}{c}{\mathrm{\&delta;F}}_x\left(t\right)={\mathrm{\&lambda;}}_{x0}+{\mathrm{\&lambda;}}_{x1}(t/10)+{\mathrm{\&lambda;}}_{x2}{(t/10)}^2+{\mathrm{\&lambda;}}_{x3}{(t/10)}^3+{\mathrm{\&lambda;}}_{x4}{(t/10)}^4+{\mathrm{\&lambda;}}_{x5}{(t/10)}^5\\ {\mathrm{\&delta;F}}_y\left(t\right)={\mathrm{\&lambda;}}_{y0}+{\mathrm{\&lambda;}}_{\mathrm{y1}}(t/10)+{\mathrm{\&lambda;}}_{y2}{(t/10)}^2+{\mathrm{\&lambda;}}_{y3}{(t/10)}^3+{\mathrm{\&lambda;}}_{y4}{(t/10)}^4+{\mathrm{\&lambda;}}_{y5}{(t/10)}^5\\ {\mathrm{\&delta;F}}_z\left(t\right)={\mathrm{\&lambda;}}_{z0}+{\mathrm{\&lambda;}}_{\mathrm{z1}}(t/10)+{\mathrm{\&lambda;}}_{z2}{(t/10)}^2+{\mathrm{\&lambda;}}_{z3}{(t/10)}^3+{\mathrm{\&lambda;}}_{z4}{(t/10)}^4+{\mathrm{\&lambda;}}_{z5}{(t/10)}^5\end{array}---\left(14\right);$

$\begin{array}{c}{\mathrm{\&delta;v}}_{ax}\left(t_k\right)=\frac{1}{n_b}{\mathrm{\&lambda;}}_{x0}\sin \left(n_b{t}_k\right)+\frac{1}{{n_b}^2}{\mathrm{\&lambda;}}_{x1}(1-\cos (n_b{t}_k)+\frac{2}{{n_b}^3}{\mathrm{\&lambda;}}_{x2}(n_b{t}_k-\sin \left(n_b{t}_k\right))+\\ \frac{3}{{n_b}^4}{\mathrm{\&lambda;}}_{x3}(-2+n_b^2{t}_k^2+2\cos (n_b{t}_k\left)\right)+\frac{4}{{n_b}^5}{\mathrm{\&lambda;}}_{x4}(-6n_b{t}_k+n_b^3{t}_k^3+6\sin (n_b{t}_k\left)\right)+\\ \frac{5}{{n_b}^6}{\mathrm{\&lambda;}}_{x5}(24-12n_b^2{t}_k^2+n_b^4{t}_k^4-24\cos (n_b{t}_k\left)\right)+\frac{1}{90720}\&CenterDot;\{n_b^2{t}_k^3\&CenterDot;\&lsqb;8064{\mathrm{\&lambda;}}_{y0}+\\ t_k\&CenterDot;(1848{\mathrm{\&lambda;}}_{\mathrm{y1}}+(696{\mathrm{\&lambda;}}_{y2}+333{\mathrm{\&lambda;}}_{\mathrm{y3}}{t}_k+184333{\mathrm{\&lambda;}}_{\mathrm{y4}}{t}_k^2+112{\mathrm{\&lambda;}}_{y5}{t}_k^3\left)\right)\&rsqb;\}\end{array}---\left(15\right);$

$\begin{array}{c}{\mathrm{\&delta;v}}_{ay}\left(t_k\right)=\frac{1}{\sqrt{2}{n}_b}{\mathrm{\&lambda;}}_{y0}\sinh \left(\sqrt{2}{n}_b{t}_k\right)+\frac{1}{2{n_b}^2}{\mathrm{\&lambda;}}_{y1}\left(\cosh \right(\sqrt{2}{n}_b{t}_k)-1)+\frac{1}{2{n_b}^3}{\mathrm{\&lambda;}}_{y2}(-2n_b{t}_k+\\ \sqrt{2}\sinh \left(\sqrt{2}{n}_b{t}_k\right))+\frac{3}{2{n_b}^4}{\mathrm{\&lambda;}}_{y3}(-1-n_b^2{t}_k^2+\cosh (\sqrt{2}{n}_b{t}_k\left)\right)+\frac{1}{{n_b}^5}{\mathrm{\&lambda;}}_{\mathrm{y4}}(-2n_b{t}_k\mathrm{(3+}\\ n_b^2{t}_k^2)+3\sqrt{2}\sinh (\sqrt{2}{n}_b{t}_k\left)\right)-\frac{5}{{n_b}^6}{\mathrm{\&lambda;}}_{y5}(6+6n_b^2{t}_k^2+n_b^4{t}_k^4-6\cosh (\sqrt{2}{n}_b{t}_k\left)\right)+\\ \frac{1}{90720}\&CenterDot;\{n_b^2{t}_k^3\&CenterDot;\&lsqb;8064{\mathrm{\&lambda;}}_{x0}+t_k\&CenterDot;(1848{\mathrm{\&lambda;}}_{\mathrm{x1}}+(696{\mathrm{\&lambda;}}_{x2}+333{\mathrm{\&lambda;}}_{\mathrm{x3}}{t}_k+184333{\mathrm{\&lambda;}}_{\mathrm{x4}}{t}_k^2+\\ 112{\mathrm{\&lambda;}}_{x5}{t}_k^3\left)\right)\&rsqb;\}\end{array}---\left(16\right);$

$\begin{array}{c}{\mathrm{\&delta;v}}_{az}\left(t_k\right)=\frac{1}{n_b}{\mathrm{\&lambda;}}_{z0}\sin \left(n_b{t}_k\right)+\frac{1}{{n_b}^2}{\mathrm{\&lambda;}}_{z1}(1-\cos (n_b{t}_k)+\frac{2}{{n_b}^3}{\mathrm{\&lambda;}}_{z3}(n_b{t}_k-\sin \left(n_b{t}_k\right))+\\ \frac{3}{{n_b}^4}{\mathrm{\&lambda;}}_{z3}(-2+n_b^2{t}_k^2+2\cos (n_b{t}_k\left)\right)+\frac{4}{{n_b}^5}{\mathrm{\&lambda;}}_{\mathrm{z4}}(-6n_b{t}_k+n_b^3{t}_k^3+6\sin (n_b{t}_k\left)\right)+\\ \frac{5}{{n_b}^6}{\mathrm{\&lambda;}}_{z5}(24-12n_b^2{t}_k^2+n_b^4{t}_k^4-24\cos (n_b{t}_k\left)\right)\end{array}---\left(17\right);$

$\begin{array}{c}{\mathrm{\&delta;x}}_a\left(t_k\right)=\frac{1}{{n_b}^7}\{n_b\&lsqb;120{\mathrm{\&lambda;}}_{x5}{t}_k+{\mathrm{\&lambda;}}_{x4}(24-12n_b^2{t}_k^2+n_b^4{t}_k^4)+n_b^2(-2{\mathrm{\&lambda;}}_{x2}-6{\mathrm{\&lambda;}}_{x3}{t}_k+n_b^2({\mathrm{\&lambda;}}_{x0}+\\ {\mathrm{\&lambda;}}_{x1}{t}_k\left)\right)+n_b^2{t}_k^2(-20{\mathrm{\&lambda;}}_{x5}{t}_k+n_b^2({\mathrm{\&lambda;}}_{x2}+{\mathrm{\&lambda;}}_{x3}{t}_k+{\mathrm{\&lambda;}}_{x5}{t}_k^3))\&rsqb;-n_b(24{\mathrm{\&lambda;}}_{x4}-2{\mathrm{\&lambda;}}_{x2}{n}_b^2+\\ {\mathrm{\&lambda;}}_{x0}{n}_b^4)\cos \left(n_b{t}_k\right)-(120{\mathrm{\&lambda;}}_{x5}-6{\mathrm{\&lambda;}}_{x3}{n}_b^2+{\mathrm{\&lambda;}}_{x1}{n}_b^4)\&CenterDot;\sin \left(n_b{t}_k\right)\}\end{array}---\left(18\right);$

$\begin{array}{c}{\mathrm{\&delta;y}}_a\left(t_k\right)=\frac{1}{2{n_b}^2}{\mathrm{\&lambda;}}_{y0}\left(\cosh \right(\sqrt{2}{n}_b{t}_k)-1)+\frac{1}{4{n_b}^3}{\mathrm{\&lambda;}}_{y1}(-2n_b{t}_k+\sqrt{2}\sinh (\sqrt{2}{n}_b{t}_k\left)\right)-\frac{1}{2{n_b}^4}{\mathrm{\&lambda;}}_{y2}(-\\ 1-n_b^2{t}_k^2+\cosh \left(\sqrt{2}{n}_b{t}_k\right))-\frac{1}{4{n_b}^5}{\mathrm{\&lambda;}}_{y3}(6n_b{t}_k+2n_b^3{t}_k^3-3\sqrt{2}\sinh (\sqrt{2}{n}_b{t}_k\left)\right)-\\ \frac{1}{2{n_b}^6}{\mathrm{\&lambda;}}_{y4}(6+6n_b^2{t}_k^2+n_b^4{t}_k^4-6\cosh (\sqrt{2}{n}_b{t}_k\left)\right)-\frac{1}{2{n_b}^7}{\mathrm{\&lambda;}}_{y5}(30n_b{t}_k+10n_b^3{t}_k^3+\\ n_b^5{t}_k^5-15\sqrt{2}\sinh \left(\sqrt{2}{n}_b{t}_k\right))\end{array}---\left(19\right);$

$\begin{array}{c}{\mathrm{\&delta;z}}_a\left(t_k\right)=\frac{1}{{n_b}^7}\{n_b\&lsqb;120{\mathrm{\&lambda;}}_{z5}{t}_k+{\mathrm{\&lambda;}}_{z4}(24-12n_b^2{t}_k^2+n_b^4{t}_k^4)+n_b^2(-2{\mathrm{\&lambda;}}_{z2}-6{\mathrm{\&lambda;}}_{z3}{t}_k+\\ n_b^2\left({\mathrm{\&lambda;}}_{z0}+{\mathrm{\&lambda;}}_{z1}{t}_k\right))+n_b^2{t}_k^2(-20{\mathrm{\&lambda;}}_{z5}{t}_k+n_b^2({\mathrm{\&lambda;}}_{z2}+{\mathrm{\&lambda;}}_{z3}{t}_k+{\mathrm{\&lambda;}}_{z5}{t}_k^3\left)\right)\&rsqb;-\\ n_b(24{\mathrm{\&lambda;}}_{z4}-2{\mathrm{\&lambda;}}_{z2}{n}_b^2+{\mathrm{\&lambda;}}_{z0}{n}_b^4)\cos \left(n_b{t}_k\right)-(120{\mathrm{\&lambda;}}_{z5}-6{\mathrm{\&lambda;}}_{z3}{n}_b^2+{\mathrm{\&lambda;}}_{z1}{n}_b^4)\&CenterDot;\sin \left(n_b{t}_k\right)\}\end{array}---\left(20\right);$

第四步，视加速度偏差量求解，具体为：根据视角加速度偏差的表达式(21)得到助推段弹道每一点视加速度偏差与该点的速度、位置矢量偏差之间的表达式(28)；

$\mathrm{\&delta;}\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{W}}=\mathrm{\&delta;}\stackrel{\&OverBar;}{R}+\mathrm{\&delta;}\stackrel{\&OverBar;}{T}---\left(21\right);$

其中和各自的表达式详见(22)，详情如下：

$\stackrel{\&OverBar;}{R}=\frac{1}{m}\left[\begin{array}{c}-C_x\frac{1}{2}\mathrm{\&rho;}{v}^2{S}_m\\ C_y\frac{1}{2}{\mathrm{\&rho;v}}^2{S}_m\\ C_z\frac{1}{2}{\mathrm{\&rho;v}}^2{S}_m\end{array}\right]---\left(22\right);$

$\stackrel{\&OverBar;}{T}=\frac{1}{m}\left[\begin{array}{c}-\stackrel{\&CenterDot;}{m}{u}_e+S_e(p_e-p_H)\\ 0\\ 0\end{array}\right]$

记M_v为对速度矢量的偏导数，M_r为对位置矢量的偏导数，N_r为对位置矢量的偏导数；

$\mathrm{\&delta;}\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{W}}=\left[\begin{array}{cc}{M}_v& M_r+N_r\end{array}\right]\left(\begin{array}{c}\mathrm{\&delta;}\stackrel{\&OverBar;}{V}\\ \mathrm{\&delta;}\stackrel{\&OverBar;}{\mathrm{\&rho;}}\end{array}\right)---\left(28\right);$

第五步，解析模型迭代修正，具体为：通过不断迭代修正扰动量使得位置、速度状态矢量偏差值逐渐逼近真实值。

3.根据权利要求2所述的扰动引力作用下弹道助推段状态偏差解析计算方法，其特征在于：所述第二步中状态偏差模型建立具体过程如下：

标准弹道计算模型表示为表达式(1)，详情如下：

$\left\{\begin{array}{c}{\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{V}}}^{*}={\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{W}}}^{*}+g^{*}\left({\stackrel{\&OverBar;}{\mathrm{\&rho;}}}^{*}\right)\\ {\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{\mathrm{\&rho;}}}}^{*}={\stackrel{\&OverBar;}{V}}^{*}\end{array}\right.---\left(1\right);$

扰动引力对导弹运动的影响时实际飞行弹道与标准弹道的等时变分表示为表达式(2)，详情如下：

$\left\{\begin{array}{c}\mathrm{\&delta;}\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{V}}=\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{W}}-{\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{W}}}^{*}+\stackrel{\&OverBar;}{g}\left(\stackrel{\&OverBar;}{\mathrm{\&rho;}}\right)-{\stackrel{\&OverBar;}{g}}^{*}\left({\stackrel{\&OverBar;}{\mathrm{\&rho;}}}^{*}\right)\\ \mathrm{\&delta;}\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{\mathrm{\&rho;}}}=\mathrm{\&delta;}\stackrel{\&OverBar;}{V}\end{array}\right.---\left(2\right);$

将表达式(2)改写为矩阵形式，同时经过小偏差处理并略去高阶小量后可得表达式(3)，详情如下：

$\left(\begin{array}{c}\mathrm{\&delta;}\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{V}}\\ \mathrm{\&delta;}\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{\mathrm{\&rho;}}}\end{array}\right)=T\left(\begin{array}{c}\mathrm{\&delta;}\stackrel{\&OverBar;}{V}\\ \mathrm{\&delta;}\stackrel{\&OverBar;}{\mathrm{\&rho;}}\end{array}\right)+\left(\begin{array}{c}\mathrm{\&delta;}\stackrel{\&OverBar;}{F}\\ 0\end{array}\right)---\left(3\right);$

其中，为扰动量，且有 为扰动引力，为视加速度偏差；式中T的表达式详见表达式(4)：

表达式(4)中A₁、A₂、A₃、A₄、A₅以及A₆六者各自的表达式详见(5)：

$\{\begin{array}{c}{A}_1=-n_b^2(1-3\frac{x^2}{r^2})\\ A_2=3n_b^2\frac{x}{r}\frac{R_0+y}{r}\\ A_3=3n_b^2\frac{x}{r}\&CenterDot;\frac{z}{r}\\ A_4=-n_b^2(1-3\frac{{(R_0+y)}^2}{r^2})\\ A_5=3n_b^2\frac{z}{r}\&CenterDot;\frac{R_0+y}{r}\\ A_6=-n_b^2(1-3\frac{z^2}{r^2})\end{array}---\left(5\right);$

经过合理简化，采用伴随方程求得表达式(3)的状态转移矩阵解析解为表达式(7)，详情如下：

$\begin{array}{c}\mathrm{\&Phi;}(t_k,t)=\mathrm{\&Phi;}\left(\mathrm{\&tau;}\right)\\ =\left[\begin{array}{cccccc}{\mathrm{cosn}}_b\mathrm{\&tau;}& {\mathrm{\&Phi;}}_{12}\left(\mathrm{\&tau;}\right)& 0& -n_b{\mathrm{sinn}}_b\mathrm{\&tau;}& 0& 0\\ {\mathrm{\&Phi;}}_{21}\left(\mathrm{\&tau;}\right)& \cosh \sqrt{2}{n}_b\mathrm{\&tau;}& 0& 0& \sqrt{2}{n}_b\sinh \sqrt{2}{n}_b\mathrm{\&tau;}& 0\\ 0& 0& {\mathrm{cosn}}_b\mathrm{\&tau;}& 0& 0& -n_b{\mathrm{sinn}}_b\mathrm{\&tau;}\\ \frac{{\mathrm{sinn}}_b\mathrm{\&tau;}}{n_b}& 0& 0& {\mathrm{cosn}}_b\mathrm{\&tau;}& 0& 0\\ 0& \frac{\sin \sqrt{2}{n}_b\mathrm{\&tau;}}{\sqrt{2}{n}_b}& 0& 0& \cosh \sqrt{2}{n}_b\mathrm{\&tau;}& 0\\ 0& 0& \frac{{\mathrm{sinn}}_b\mathrm{\&tau;}}{n_b}& 0& 0& {\mathrm{cosn}}_b\mathrm{\&tau;}\end{array}\right]\end{array}---\left(7\right);$

其中Φ₁₂(τ)详见表达式(8)：

${\mathrm{\&Phi;}}_{12}\left(\mathrm{\&tau;}\right)={\mathrm{\&Phi;}}_{21}\left(\mathrm{\&tau;}\right)=\frac{{\left(n_b\mathrm{\&tau;}\right)}^2}{3!}(1+\frac{2}{3}\frac{\mathrm{\&tau;}}{t_k}+\frac{1}{6}\frac{{\mathrm{\&tau;}}^2}{t_k^2})---\left(8\right);$

根据表达式(7)即可推得助推段扰动引力引起的状态偏差半解析表达式(9)，如下：

${\mathrm{\&delta;v}}_{ax}\left(t_k\right)={\&Integral;}_0^{t_k}{\mathrm{\&Phi;}}_{11}(t_k-\mathrm{\&tau;}){\mathrm{\&delta;F}}_x\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;}+{\&Integral;}_0^{t_k}{\mathrm{\&Phi;}}_{12}(t_k-\mathrm{\&tau;}){\mathrm{\&delta;F}}_y\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;}$

${\mathrm{\&delta;v}}_{ay}\left(t_k\right)={\&Integral;}_0^{t_k}{\mathrm{\&Phi;}}_{22}(t_k-\mathrm{\&tau;}){\mathrm{\&delta;F}}_y\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;}+{\&Integral;}_0^{t_k}{\mathrm{\&Phi;}}_{21}(t_k-\mathrm{\&tau;}){\mathrm{\&delta;F}}_x\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;}$

${\mathrm{\&delta;v}}_{az}\left(t_k\right)={\&Integral;}_0^{t_k}{\mathrm{\&Phi;}}_{33}(t_k-\mathrm{\&tau;}){\mathrm{\&delta;F}}_z\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;}$

${\mathrm{\&delta;x}}_a\left(t_k\right)={\&Integral;}_0^{t_k}{\mathrm{\&Phi;}}_{41}(t_k-\mathrm{\&tau;}){\mathrm{\&delta;F}}_x\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;}---\left(9\right).$

${\mathrm{\&delta;y}}_a\left(t_k\right)={\&Integral;}_0^{t_k}{\mathrm{\&Phi;}}_{52}(t_k-\mathrm{\&tau;}){\mathrm{\&delta;F}}_y\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;}$

${\mathrm{\&delta;z}}_a\left(t_k\right)={\&Integral;}_0^{t_k}{\mathrm{\&Phi;}}_{63}(t_k-\mathrm{\&tau;}){\mathrm{\&delta;F}}_z\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;}$

4.根据权利要求3所述的扰动引力作用下弹道助推段状态偏差解析计算方法，其特征在于：所述合理简化的过程具体为：

由表达式(5)可知，A₁、A₂、A₃、A₄、A₅以及A₆六者均为小量，其中A₁、A₄以及A₆较大，其数量级与相当，而A₂比小一个数量级，A₃和A₅则比小两三个数量级，即有：

$\left\{\begin{array}{c}\frac{R_0+y}{r}\&ap;1\\ \frac{z^2}{r^3}<<1\\ \frac{x^2}{r^2}<<1\end{array}\right.;$

因此，将表达式(5)改写为表达式(6)的形式，详情如下：

$\left\{\begin{array}{c}{A}_1=A_6=-n_b^2\\ A_4=2n_b^2\\ A_2=\frac{3x}{r}{n}_b^2\\ A_3=A_5=0\end{array}\right.---\left(6\right);$

其中，n_b为常数。

5.根据权利要求3所述的扰动引力作用下弹道助推段状态偏差解析计算方法，其特征在于：所述第三步中助推段扰动量拟合具体为：

采用最小二乘法将助推段扰动量拟合为关于时间的多项式函数，且发射惯性系中扰动量三分量分别进行拟合，具体是：设观测方程为表达式(10)：

Z_n×1＝H_n×(m+1)λ_(m+1)×1+ε(10)；

其中，Z_n×1为观测向量，n为拟合点的个数；H_n×(m+1)为系数矩阵，m代表了多项式拟合次数；λ_(m+1)×1为待定参数向量；ε为随机误差向量；

拟合多项式采用经验公式(11)则可得系数矩阵和扰动量拟合系数为表述式(12)和(13)，详情如下：

$\mathrm{\&delta;}F\left(t\right)=\underset{i=0}{\overset{n}{\&Sigma;}}{\mathrm{\&lambda;}}_i{\left(\frac{t}{10}\right)}^i---\left(11\right);$

λ＝(H^TH)^-1H^TZ(13)；

当拟合多项式取5阶时，可得扰动量拟合表达式(14)，如下：

$\left\{\begin{array}{c}{\mathrm{\&delta;F}}_x\left(t\right)={\mathrm{\&lambda;}}_{x0}+{\mathrm{\&lambda;}}_{x1}(t/10)+{\mathrm{\&lambda;}}_{x2}{(t/10)}^2+{\mathrm{\&lambda;}}_{x3}{(t/10)}^3+{\mathrm{\&lambda;}}_{x4}{(t/10)}^4+{\mathrm{\&lambda;}}_{x5}{(t/10)}^5\\ {\mathrm{\&delta;F}}_y\left(t\right)={\mathrm{\&lambda;}}_{y0}+{\mathrm{\&lambda;}}_{y1}(t/10)+{\mathrm{\&lambda;}}_{y2}{(t/10)}^2+{\mathrm{\&lambda;}}_{y3}{(t/10)}^3+{\mathrm{\&lambda;}}_{y4}{(t/10)}^4+{\mathrm{\&lambda;}}_{y5}{(t/10)}^5\\ {\mathrm{\&delta;F}}_z\left(t\right)={\mathrm{\&lambda;}}_{z0}+{\mathrm{\&lambda;}}_{z1}(t/10)+{\mathrm{\&lambda;}}_{z2}{(t/10)}^2+{\mathrm{\&lambda;}}_{z3}{(t/10)}^3+{\mathrm{\&lambda;}}_{z4}{(t/10)}^4+{\mathrm{\&lambda;}}_{z5}{(t/10)}^5\end{array}\right.---\left(14\right);$

将表达式(14)带入助推段扰动引力引起的状态偏差半解析表达式(9)导出助推段扰动引力引起的状态偏差完整解析表达式(15)、(16)、(17)、(18)、(19)以及(20)，详情如下：

$\begin{array}{c}{\mathrm{\&delta;v}}_{ax}\left(t_k\right)=\frac{1}{n_b}{\mathrm{\&lambda;}}_{x0}\sin \left(n_b{t}_k\right)+\frac{1}{{n_b}^2}{\mathrm{\&lambda;}}_{x1}(1-\cos (n_b{t}_k)+\frac{2}{{n_b}^3}{\mathrm{\&lambda;}}_{x2}(n_b{t}_k-\sin \left(n_b{t}_k\right))+\\ \frac{3}{{n_b}^4}{\mathrm{\&lambda;}}_{x3}(-2+n_b^2{t}_k^2+2\cos (n_b{t}_k\left)\right)+\frac{4}{{n_b}^5}{\mathrm{\&lambda;}}_{x4}(-6n_b{t}_k+n_b^3{t}_k^3+6\sin (n_b{t}_k\left)\right)+\\ \frac{5}{{n_b}^6}{\mathrm{\&lambda;}}_{x5}(24-12n_b^2{t}_k^2+n_b^4{t}_k^4-24\cos (n_b{t}_k\left)\right)+\frac{1}{90720}\&CenterDot;\{n_b^2{t}_k^3\&CenterDot;\&lsqb;8064{\mathrm{\&lambda;}}_{y0}+\\ t_k\&CenterDot;(1848{\mathrm{\&lambda;}}_{\mathrm{y1}}+(696{\mathrm{\&lambda;}}_{y2}+333{\mathrm{\&lambda;}}_{\mathrm{y3}}{t}_k+184333{\mathrm{\&lambda;}}_{\mathrm{y4}}{t}_k^2+112{\mathrm{\&lambda;}}_{y5}{t}_k^3\left)\right)\&rsqb;\}\end{array}---\left(15\right);$

$\begin{array}{c}{\mathrm{\&delta;v}}_{ay}\left(t_k\right)=\frac{1}{\sqrt{2}{n}_b}{\mathrm{\&lambda;}}_{y0}\sinh \left(\sqrt{2}{n}_b{t}_k\right)+\frac{1}{2{n_b}^2}{\mathrm{\&lambda;}}_{y1}\left(\cosh \right(\sqrt{2}{n}_b{t}_k)-1)+\frac{1}{2{n_b}^3}{\mathrm{\&lambda;}}_{y2}(-2n_b{t}_k+\\ \sqrt{2}\sinh \left(\sqrt{2}{n}_b{t}_k\right))+\frac{3}{2{n_b}^4}{\mathrm{\&lambda;}}_{y3}(-1-n_b^2{t}_k^2+\cosh (\sqrt{2}{n}_b{t}_k\left)\right)+\frac{1}{{n_b}^5}{\mathrm{\&lambda;}}_{\mathrm{y4}}(-2n_b{t}_k\mathrm{(3+}\\ n_b^2{t}_k^2)+3\sqrt{2}\sinh (\sqrt{2}{n}_b{t}_k\left)\right)-\frac{5}{{n_b}^6}{\mathrm{\&lambda;}}_{y5}(6+6n_b^2{t}_k^2+n_b^4{t}_k^4-6\cosh (\sqrt{2}{n}_b{t}_k\left)\right)+\\ \frac{1}{90720}\&CenterDot;\{n_b^2{t}_k^3\&CenterDot;\&lsqb;8064{\mathrm{\&lambda;}}_{x0}+t_k\&CenterDot;(1848{\mathrm{\&lambda;}}_{\mathrm{x1}}+(696{\mathrm{\&lambda;}}_{x2}+333{\mathrm{\&lambda;}}_{\mathrm{x3}}{t}_k+184333{\mathrm{\&lambda;}}_{\mathrm{x4}}{t}_k^2+\\ 112{\mathrm{\&lambda;}}_{x5}{t}_k^3\left)\right)\&rsqb;\}\end{array}---\left(16\right);$

$\begin{array}{c}{\mathrm{\&delta;v}}_{az}\left(t_k\right)=\frac{1}{n_b}{\mathrm{\&lambda;}}_{z0}\sin \left(n_b{t}_k\right)+\frac{1}{{n_b}^2}{\mathrm{\&lambda;}}_{z1}(1-\cos (n_b{t}_k)+\frac{2}{{n_b}^3}{\mathrm{\&lambda;}}_{z3}(n_b{t}_k-\sin \left(n_b{t}_k\right))+\\ \frac{3}{{n_b}^4}{\mathrm{\&lambda;}}_{z3}(-2+n_b^2{t}_k^2+2\cos (n_b{t}_k\left)\right)+\frac{4}{{n_b}^5}{\mathrm{\&lambda;}}_{\mathrm{z4}}(-6n_b{t}_k+n_b^3{t}_k^3+6\sin (n_b{t}_k\left)\right)+\\ \frac{5}{{n_b}^6}{\mathrm{\&lambda;}}_{z5}(24-12n_b^2{t}_k^2+n_b^4{t}_k^4-24\cos (n_b{t}_k\left)\right)\end{array}---\left(17\right);$

$\begin{array}{c}{\mathrm{\&delta;x}}_a\left(t_k\right)=\frac{1}{{n_b}^7}\{n_b\&lsqb;120{\mathrm{\&lambda;}}_{x5}{t}_k+{\mathrm{\&lambda;}}_{x4}(24-12n_b^2{t}_k^2+n_b^4{t}_k^4)+n_b^2(-2{\mathrm{\&lambda;}}_{x2}-6{\mathrm{\&lambda;}}_{x3}{t}_k+n_b^2({\mathrm{\&lambda;}}_{x0}+\\ {\mathrm{\&lambda;}}_{x1}{t}_k\left)\right)+n_b^2{t}_k^2(-20{\mathrm{\&lambda;}}_{x5}{t}_k+n_b^2({\mathrm{\&lambda;}}_{x2}+{\mathrm{\&lambda;}}_{x3}{t}_k+{\mathrm{\&lambda;}}_{x5}{t}_k^3\left)\right)\&rsqb;-n_b(24{\mathrm{\&lambda;}}_{x4}-2{\mathrm{\&lambda;}}_{x2}{n}_b^2+\\ {\mathrm{\&lambda;}}_{x0}{n}_b^4)\cos \left(n_b{t}_k\right)-(120{\mathrm{\&lambda;}}_{x5}-6{\mathrm{\&lambda;}}_{x3}{n}_b^2+{\mathrm{\&lambda;}}_{x1}{n}_b^4)\&CenterDot;\sin \left(n_b{t}_k\right)\}\end{array}---\left(18\right);$

$\begin{array}{c}{\mathrm{\&delta;y}}_a\left(t_k\right)=\frac{1}{2{n_b}^2}{\mathrm{\&lambda;}}_{y0}\left(\cosh \right(\sqrt{2}{n}_b{t}_k)-1)+\frac{1}{4{n_b}^3}{\mathrm{\&lambda;}}_{y1}(-2n_b{t}_k+\sqrt{2}\sinh (\sqrt{2}{n}_b{t}_k\left)\right)-\frac{1}{2{n_b}^4}{\mathrm{\&lambda;}}_{y2}(-\\ 1-n_b^2{t}_k^2+\cosh \left(\sqrt{2}{n}_b{t}_k\right))-\frac{1}{4{n_b}^5}{\mathrm{\&lambda;}}_{y3}(6n_b{t}_k+2n_b^3{t}_k^3-3\sqrt{2}\sinh (\sqrt{2}{n}_b{t}_k\left)\right)-\\ \frac{1}{2{n_b}^6}{\mathrm{\&lambda;}}_{y4}(6+6n_b^2{t}_k^2+n_b^4{t}_k^4-6\cosh (\sqrt{2}{n}_b{t}_k\left)\right)-\frac{1}{2{n_b}^7}{\mathrm{\&lambda;}}_{y5}(30n_b{t}_k+10n_b^3{t}_k^3+\\ n_b^5{t}_k^5-15\sqrt{2}\sinh \left(\sqrt{2}{n}_b{t}_k\right))\end{array}---\left(19\right);$

$\begin{array}{c}{\mathrm{\&delta;z}}_a\left(t_k\right)=\frac{1}{{n_b}^7}\{n_b\&lsqb;120{\mathrm{\&lambda;}}_{z5}{t}_k+{\mathrm{\&lambda;}}_{z4}(24-12n_b^2{t}_k^2+n_b^4{t}_k^4)+n_b^2(-2{\mathrm{\&lambda;}}_{z2}-6{\mathrm{\&lambda;}}_{z3}{t}_k+\\ n_b^2({\mathrm{\&lambda;}}_{z0}+{\mathrm{\&lambda;}}_{z1}{t}_k))+n_b^2{t}_k^2(-20{\mathrm{\&lambda;}}_{z5}{t}_k+n_b^2({\mathrm{\&lambda;}}_{z2}+{\mathrm{\&lambda;}}_{z3}{t}_k+{\mathrm{\&lambda;}}_{z5}{t}_k^3))\&rsqb;-\\ n_b(24{\mathrm{\&lambda;}}_{z4}-2{\mathrm{\&lambda;}}_{z2}{n}_b^2+{\mathrm{\&lambda;}}_{z0}{n}_b^4)\cos \left(n_b{t}_k\right)-(120{\mathrm{\&lambda;}}_{z5}-6{\mathrm{\&lambda;}}_{z3}{n}_b^2+{\mathrm{\&lambda;}}_{z1}{n}_b^4)\&CenterDot;\sin \left(n_b{t}_k\right)\}\end{array}---\left(20\right);$

6.根据权利要求5所述的扰动引力作用下弹道助推段状态偏差解析计算方法，其特征在于：所述观测向量为基于弹道助推段n个离散时间点对应的位置矢量求解得到的n组发射坐标系中扰动引力值，所述离散时间点按照前密后松的原则进行选取。

7.根据权利要求5所述的扰动引力作用下弹道助推段状态偏差解析计算方法，其特征在于：所述第四步中视加速度偏差量求解具体为：

视加速度由气动力和发动机推力两部分组成，故而视角加速度偏差可表示为表达式(21)，详情如下： $\mathrm{\&delta;}\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{W}}=\mathrm{\&delta;}\stackrel{\&OverBar;}{R}+\mathrm{\&delta;}\stackrel{\&OverBar;}{T}---\left(21\right);$

其中和各自的表达式详见(22)，详情如下：

$\stackrel{\&OverBar;}{R}=\frac{1}{m}\left[\begin{array}{c}-C_x\frac{1}{2}\mathrm{\&rho;}r{v}^2{S}_m\\ C_y\frac{1}{2}{\mathrm{\&rho;v}}^2{S}_m\\ C_z\frac{1}{2}{\mathrm{\&rho;v}}^2{S}_m\end{array}\right]---\left(22\right);$

$\stackrel{\&OverBar;}{T}=\frac{1}{m}\left[\begin{array}{c}-\stackrel{\&CenterDot;}{m}{u}_e+S_e(p_e-p_H)\\ 0\\ 0\end{array}\right]$

记M_v为对速度矢量的偏导数，M_r为对位置矢量的偏导数，具体表示如表达式(23)，详情如下：

$M_v=\frac{\&part;\stackrel{\&OverBar;}{R}}{\&part;{\stackrel{\&OverBar;}{V}}^T}=\frac{S_m}{m}\left[\begin{array}{ccc}-C_x{\mathrm{\&rho;v}}_x& -C_x{\mathrm{\&rho;v}}_y& -C_x{\mathrm{\&rho;v}}_z\\ C_y^{\mathrm{\&alpha;}}\mathrm{\&rho;}({\mathrm{\&alpha;v}}_x+\frac{v^2{V}_y}{2(V_x^2+V_y^2)})& C_y^{\mathrm{\&alpha;}}\mathrm{\&rho;}({\mathrm{\&alpha;v}}_y-\frac{v^2{V}_x}{2(V_x^2+V_y^2)})& C_y^{\mathrm{\&alpha;}}{\mathrm{\&rho;\&alpha;v}}_z\\ C_z^{\mathrm{\&beta;}}\mathrm{\&rho;}({\mathrm{\&beta;v}}_x-\frac{V_x{V}_y}{2\sqrt{V_x^2+V_y^2}})& C_z^{\mathrm{\&beta;}}\mathrm{\&rho;}({\mathrm{\&beta;v}}_y+\frac{V_y{V}_z}{2\sqrt{V_x^2+V_y^2}})& C_z^{\mathrm{\&beta;}}\mathrm{\&rho;}({\mathrm{\&beta;v}}_z+\frac{v^2}{2\sqrt{V_x^2+V_y^2}})\end{array}\right]---\left(23\right);$

$M_r=\frac{\&part;\stackrel{\&OverBar;}{R}}{\&part;{\stackrel{\&OverBar;}{\mathrm{\&rho;}}}^T}=\frac{S_m}{m}\left[\begin{array}{ccc}0& -C_x\frac{1}{2}{v}^2\frac{\&part;\mathrm{\&rho;}}{\&part;y}& 0\\ 0& C_y^{\mathrm{\&alpha;}}\mathrm{\&alpha;}\frac{1}{2}{v}^2\frac{\&part;\mathrm{\&rho;}}{\&part;y}& 0\\ 0& C_z^{\mathrm{\&beta;}}\mathrm{\&beta;}\frac{1}{2}{v}^2\frac{\&part;\mathrm{\&rho;}}{\&part;y}& 0\end{array}\right]$

其中，为大气密度对高度的偏导数；

记N_r为对位置矢量的偏导数，具体详见表达式(24)：

$N_{\mathrm{\&gamma;}}=\frac{\&part;\stackrel{\&OverBar;}{T}}{\&part;{\stackrel{\&OverBar;}{\mathrm{\&rho;}}}^T}=\frac{S_e}{m}\left[\begin{array}{ccc}0& -\frac{\&part;p_H}{\&part;y}& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]---\left(24\right);$

其中，为大气压强对高度的偏导数；

弹体系到发射系、速度系到发射系以及发射惯性系到发射系的转换矩阵分别为表达式(25)、(26)以及(27)：

$G_A=\left[\begin{array}{ccc}1& {\mathrm{\&omega;}}_{ez}t& -{\mathrm{\&omega;}}_{ey}t\\ -{\mathrm{\&omega;}}_{ex}t& 1& {\mathrm{\&omega;}}_{ez}t\\ {\mathrm{\&omega;}}_{ey}t& -{\mathrm{\&omega;}}_{ex}t& 1\end{array}\right]---\left(27\right);$

其中，G_b中忽略了滚转角，G_v中忽略了倾侧角，G_A中各项近似至ω_et的一次项；

将表达式(23)、(24)代入表达式(21)中得到助推段弹道每一点视加速度偏差与该点的速度、位置矢量偏差之间的表达式(28)：

$\mathrm{\&delta;}\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{W}}=\left[\begin{array}{cc}{M}_v& M_r+N_r\end{array}\right]\left(\begin{array}{c}\mathrm{\&delta;}\stackrel{\&OverBar;}{V}\\ \mathrm{\&delta;}\stackrel{\&OverBar;}{\mathrm{\&rho;}}\end{array}\right)---\left(28\right).$

8.根据权利要求7所述的扰动引力作用下弹道助推段状态偏差解析计算方法，其特征在于：所述第五步中解析模型迭代修正的具体迭代过程如下：

迭代初值为：取扰动量即取0；求解扰动引力并保存到数组中，迭代过程中不需重新计算；同时取

步骤一：对扰动量进行拟合；

步骤二：由表达式(15)、(16)、(17)、(18)、(19)以及(20)计算并输出从开始到发动机关机整个过程中每个时间间隔的状态偏差 ${\left[\begin{array}{cc}\mathrm{\&delta;}\stackrel{\&OverBar;}{V}\left(t\right)& \mathrm{\&delta;}\stackrel{\&OverBar;}{\mathrm{\&rho;}}\left(t\right)\end{array}\right]}^T,$ 其中输出间隔与扰动量拟合取样间隔一致；

步骤三：由表达式(28)计算步骤二中每一组输出量对应的视加速度偏差

步骤四：更新每一点扰动量的值，即

步骤五：判断 $\left\{\begin{array}{c}\left|\right|\mathrm{\&delta;}\stackrel{\&OverBar;}{V}(t_k,n)-\mathrm{\&delta;}\stackrel{\&OverBar;}{V}(t_k,n-1)\left|\right|<=eps\_v\\ \left|\right|\mathrm{\&delta;}\stackrel{\&OverBar;}{\mathrm{\&rho;}}(t_k,n)-\mathrm{\&delta;}\stackrel{\&OverBar;}{\mathrm{\&rho;}}(t_k,n-1)\left|\right|<=eps\_r\end{array}\right.$ 是否成立，即相邻两次关机点状态偏差是否小于一定值，若成立，则迭代结束；若不成立，则返回步骤一。