CN105701283B - 地球非球形摄动作用下自由段弹道误差传播的分析方法 - Google Patents

Abstract

本发明提供了一种地球非球形摄动作用下自由段弹道误差传播的分析方法，针对典型的弹道导弹，首先在考虑地球非球形摄动影响的情况下，建立适用于自由段运动分析与计算的高精度简化模型，然后根据小偏差假设与线性系统理论，导出了等角/等地心距/等时摄动模型及其状态转移矩阵解析解，进一步利用纵向和侧向摄动方程的解耦，基于最小二乘拟合方法对高阶项引起的侧向偏差进行修正，并实现地球非球形摄动影响下自由段弹道特性计算的快速性与精确性的平衡。本发明方法具有计算速度快、适应区域广、修正模型精的特征，能满足快速机动发射以及弹上实时计算的要求。

Description

地球非球形摄动作用下自由段弹道误差传播的分析方法

技术领域

本发明涉及飞行动力学技术领域，具体涉及一种地球非球形摄动作用下自由段弹道误差传播的分析方法。

背景技术

引起弹道导弹落点偏差的主要因素包括制导工具误差和制导方法误差。随着惯性测量系统硬件水平的提高，可有效修正部分制导工具误差，使得弹道导弹制导工具误差逐渐降低，从而使得导方法误差的影响也就日益突出，况且地球非球形摄动正是引起制导方法误差的主要因素。对射程一万多公里的洲际弹道导弹而言，地球扁率摄动对导弹落点的影响可达几十公里，而全程扰动引力场也可达数百米甚至公里量级。

即使在少数情况下，对地球非球形摄动因素进行了初步补偿，部分减少了非球形摄动引起的命中精度偏差，但由于其地球非球形摄动因素的表征理论不完备、误差传播特性不明确、补偿方法存在各种缺陷(如补偿阶次很低，一般只补偿带谐项，如J2/J4/J6等前几项2补偿方式一般采用事先补偿，即将补偿量作为诸元装订到弹上，而并非实时计算补偿量)，致使补偿后的非球形摄动因素引起的落点偏差依然很大。

因此，为了实现弹道导弹的快速机动发射和精确打击能力，地球非球形摄动因素影响的弹道误差传播理论与快速精确补偿方法已成为亟待解决的关键问题。其中，地球非球形摄动影响的弹道误差传播理论更是快速精确补偿的前提和基础，对有效提高弹道导弹命中精度具有重要的意义。该弹道误差传播方法，旨在解决当前弹道误差传播理论领域普遍存在的三大问题：一是，建立适用于自由段弹道分析与计算的高精度简化模型；二是，提出描述自由段弹道误差传播特性的系统分析方法；三是，实现地球非球形摄动影响下自由段弹道特性计算的快速性与精确精确性的平衡。

发明内容

本发明目的在于提供一种地球非球形摄动作用下自由段弹道误差传播的分析方法，包括以下步骤：

步骤一：选定自由段射程角序列，具体是：根据一定的自由段射程角间隔△β和自由段射程角β，确定一组自由段射程角序列{β_i}；

步骤二：纵向摄动方程和侧向摄动方程的获得，具体是：建立以自由段标准射程角β为自变量的弹道导弹轨道柱坐标系简化运动模型，基于将简化的自由段运动方程，借助弹道摄动思想将此非线性方程线性化，并导出等角/等地心距/等时摄动模型及其状态转移矩阵解析解得出其状态转移矩阵解析解，包括：a、在一条标准弹道附近，将关于摄动偏差量的弹道导弹自由段非线性运动方程进行线性化，得到摄动模型下的摄动状态方程；b、根据线性系统基本理论，推导该线性时变系统状态转移矩阵的解析解；c、通过状态转移矩阵的解析解得到纵向摄动方程和侧向摄动方程；

步骤三：经过高阶项的最小二乘拟合后获得高阶摄动项的解析表达式，具体是：首先，设定观测方程和选定7阶拟合多项式作为拟合公式，其次，根据最小二乘法得到待定参数向量的线性无偏最优估计表达式；最后，获得高阶项的最小二乘估计式表达式，即为高阶摄动项的解析表达式；

步骤四：对状态参数进行高阶偏差修正，具体是：首先，取定当自由段标准射程角为β时高阶修正的过程偏差表达式；其次，对地心距偏差项进行高阶修正；最后，得到高阶修正的纵向偏差参数表达式。

以上技术方案中优选的，所述步骤一中：所述自由段射程角β为从主动段关机点起算的标准轨道射程角，所有真实轨道的自由段射程角均需投影至标准轨道面来进行确定；所述自由段射程角间隔△β为0.01度。

以上技术方案中优选的，所述步骤二中自由段运动方程的简化过程如下：

在轨道柱坐标系中，三个坐标轴的单位矢量表示为表达式1)：

其中，β为自由段标准射程角；

三个坐标轴的单位矢量的微分分别为：

以自由段起点为时间原点，任一时刻t的真实地心矢量由该时刻的标准轨道地心距r和标准侧向位移z表示如表达式3)：

则其速度表示为表达式4)：其中为r的导数；

定义沿三个坐标轴方向的标准径向速度V_r、标准周向速度V_β和标准侧向速度V_z分别为表达式5)：

则其加速度表示为表达式6)：

其中为V_z的导数；

在地球非球型摄动加速度影响的情况下，运用牛顿第二定律建立起受力状态与运动状态之间的关系如表达式7)：

其中，μ为地球引力常数；η为标准侧向偏差角，且δa_r、δa_β和δa_z分别为摄动加速度投影至轨道柱坐标系的三分量；

将标准侧向偏差角η和标准侧向速度V_z均看作小量，略去二阶及以上小量，得到以自由段标准射程角β为自变量的弹道导弹轨道柱坐标系简化运动模型，详情如表达式8)：

其中，上标“'”表示对自由段标准射程角取微分；

对于地球非球形摄动，相应的摄动加速度表示为非球形摄动引力位的梯度，再将其投影至标准轨道坐标系，得到分量形式的表达式9)：

以上技术方案中优选的，所述步骤二中的摄动状态方程具体为表达式10)：

其中，δX为偏差状态量，V为摄动项。

以上技术方案中优选的，所述步骤二中的摄动模型为等角摄动模型、等地心摄动模型以及等时摄动模型，三种摄动模型的偏差状态Y为表达式12)：

Y＝[y₁ y₂ y₃ y₄ y₅ y₆]^Τ 12)，

以标准弹道为参考弹道，对简化的弹道导弹轨道柱坐标系自由段运动模型进行线性化，得三种摄动模型下的摄动状态方程如下：

等角摄动状态方程如表达式13)：

其中，

等地心距摄动状态方程如表达式14)：

其中：

等时摄动状态方程如表达式15)：

其中：

以上技术方案中优选的，所述等角摄动状态转移矩阵的解析解的推导过程如下：

令

则等角摄动状态方程转变为表达式17)：

同时，给定初始条件

x_i(β₀)＝x_i(β₀)i＝1,2,…,6 18)，

将表达式17)中的第三式求解得到表达式19)：x₃(β)＝x₃(β₀) 19)，

将表达式17)中的第三式代入第一式中即得到表达式20)：

将表达式20)代入表达式17)中的第二式得表达式21)：

把表达式17)中的第六式代入第五式将侧向偏差参数解出如表达式22)：

将式表达式19)和表达式21)代入表达式17)中的第四式，完成纵向偏差参数的求解如表达式23)：

式中，

其中，f、r、p、h和e分别为真近点角、地心距、半通径、单位质量动量矩和偏心率；

将纵向偏差参数求解表达式23)整理为矩阵形式如表达式24)：

Φ(β,β₀)＝[Φ_ij(β,β₀)]i,j＝1,2,…,6 24)，

其中，

Φ₁₁(β,β₀)＝cos(β-β₀)；

Φ₁₂(β,β₀)＝-sin(β-β₀)；

Φ₂₁(β,β₀)＝sin(β-β₀)；

Φ₂₂(β,β₀)＝cos(β-β₀)；

Φ₄₄(β,β₀)＝1；

Φ₅₅(β,β₀)＝cos(β-β₀)；

Φ₅₆(β,β₀)＝-sin(β-β₀)；

Φ₆₅(β,β₀)＝sin(β-β₀)；

Φ₆₆(β,β₀)＝cos(β-β₀)；

其余未列出各项均为0。

以上技术方案中优选的，所述步骤二中：矩阵形式的等地心距摄动状态转移矩阵解析解表示为表达式25)：

λ_r(β,β₀)＝[λ_rij(β,β₀)]i,j＝1,2,…,6 25)，

其中：

λ_r22(β,β₀)＝1；

λ_r44(β,β₀)＝1；

其余未列出各项均为0；

矩阵形式的等时摄动状态转移矩阵解析解表示为表达式26)：

λ_t(β,β₀)＝[λ_tij(β,β₀)]i,j＝1,2,…,6 26)，

其中：

λ_t44(β,β₀)＝1；

式中，

其余未列出各项均为0。

以上技术方案中优选的，当自由段标准射程角为β时，等角摄动、等地心距摄动以及等时摄动三种摄动模型的过程偏差均表示为表达式27)：

其中，U(ξ)为等角摄动/等地心距摄动/等时过程摄动项，具体如表达式28)：

将表达式28)代入表达式27)得纵向摄动方程如表达式29)以及侧向摄动方程如表达式30)：

其中，μ为地球引力常数，p为标准椭圆弹道的半通径，r为地心距，y_ii＝1…4为纵向摄动偏差，δa_r和δa_β为摄动加速度投影在径向和周向的分量，λ为等角摄动/等地心距摄动/等时摄动状态转移矩阵解析解；

其中，μ为地球引力常数，p为标准椭圆弹道的半 通径，r为地心距，y_i i＝5,6为侧向摄动偏差，δa_r为摄动加速度投影在径向的分量，λ为等 角/等地心距/等时摄动状态转移矩阵解析解。

以上技术方案中优选的，所述步骤三中：所述观测方程为表达式31)：

z＝ha+ε 31)，

其中，z为观测向量，即以一定自由段标准射程角间隔积分解算半解析的侧向摄动方程所得到的一组侧向偏差角；h为系数矩阵，由线性拟合公式确定；a为待定参数向量；ε为随机误差向量；

所述7阶拟合多项式作为拟合公式如表达式33)：

η²＝a₀+a₁sinβ+a₂cosβ+a₃sin2β+a₄cos2β+a₅sin3β+a₆cos3β 33)，

待定参数向量a与系数矩阵h详见表达式34)和表达式37)，观测向量z与随机误差向量ε的表达式详见表达式35)和表达式36)：

a＝[a₀ a₁ a₂ a₃ a₄ a₅ a₆]^Τ 34)；

ε＝[ε₁ ε₂ … ε_n]^Τ 36)；

且

所述待定参数向量的线性无偏最优估计表达式如表达式38)：

a＝(h^Th)^-1h^Tz 38)；

所述高阶项的最小二乘估计式表达式如表达式39)：

以上技术方案中优选的，所述步骤四中对等地心距摄动的状态参数进行高阶偏差修正的具体过程如下：

当自由段标准射程角为β时，等地心距摄动高阶修正的过程偏差表示为表达式40)：

其中，W(β₀)为等地心距过程摄动项，△U_r(ξ)为高阶修正项，具体如表达式41)和表达式42)：

U_r(ξ)＝[△u_r 0 0 0 0 0]^Τ 42)；

高阶摄动项的最小二乘修正法是为了修正径向参数，因此，地心距偏差项的高阶修正为表达式43)：

将等地心距摄动状态转移矩阵解析解中λ_r21(β,ξ)的表达式代入表达式43)即得高阶修正的等地心摄动纵向偏差参数如表达式44)：

其中：

u₀(β)＝1-cosβ；

所述步骤四中高阶修正的等角摄动纵向偏差参数详见表达式45)：

其中，其中：

u₀(β)＝1-cosβ；

本发明针对典型的弹道导弹，首先在考虑地球非球形摄动影响的情况下，建立适用于自由段运动分析与计算的高精度简化模型，然后根据小偏差假设与线性系统理论，导出了等角/等地心距/等时摄动模型及其状态转移矩阵解析解，进一步利用纵向和侧向摄动方程的解耦，基于最小二乘拟合方法对高阶项引起的侧向偏差进行修正，并实现地球非球形摄动影响下自由段弹道特性计算的快速性与精确性的平衡。本发明方法具有计算速度快、适应区域广、修正模型精的特征，能满足快速机动发射以及弹上实时计算的要求，详情如下：

(1)本发明基于小偏差理论和线性系统理论构建了自由段摄动模型，模型简化合理，易于导出状态转移矩阵的解析解。

(2)本发明可结合各种地球非球形摄动模型进行误差传播分析与计算，且模型精度与复杂性可根据需要任意调整，具有良好的适应性、开放性和可扩展性。

(3)本发明考虑了自由段摄动模型中的一阶主项，并修正了线性化过程中所忽略的高阶项主项，与弹道求差法相比，修正后地球扁率J2项影响的误差传播分析误差优于20m(补偿了地球扁率J2项引起误差的99.8％)，扰动引力场影响的误差传播误差优于3m(补偿了扰动引力场引起误差的98％)。

(4)本发明具有计算速度快、适应区域广、修正模型精、弹上存储量小等的特征，状态偏差的计算具有实时性，模型计算时间主要由地球非球形摄动模型的计算时间决定，且该计算方法适应于任意弹道。

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

附图说明

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

图1为标准轨道坐标系示意图；

图2为实施例1中弹道摄动思想的框架图；

图3为实施例1中等地心距摄动偏差与等角摄动偏差之间的关系；

图4为实施例1中5项、7项以及9项拟合公式的拟合精度对比图；

图5为J2项引力为引起的位置误差示意图；

图6为J2项引力为引起的位置误差进行高阶修正后的效果示意图；

图7为36阶扰动引力场引起的位置误差示意图；

图8为36阶扰动引力场引起的位置误差进行高阶修正后的效果示意图；

其中，a、真实轨道，b、标准轨道面，c、实际弹道，d、标准椭圆弹道，e、等角线，f、等地心距线。

具体实施方式

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

实施例1：

一种地球非球形摄动作用下自由段弹道误差传播的分析方法，详情如下：

分别取摄动项为J2项引力加速度和36阶扰动引力加速度(详见图5、图6、图7和图8)，以地心距为6465411m当作终端条件，对发射方位角进行全方位循环，以验证本发明分析计算的精度。

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

引力模型为：EGM2008模型(EGM指地球重力场)。

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

三个摄动模型：等地心距摄动模型、等角摄动模型以及等时摄动模型。

具体包含以下步骤：

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

表1仿真初始参数

分析方法具体包含以下步骤：

第一步，选定自由段射程角序列，具体是：根据一定的自由段射程角间隔△β和自由段射程角β，确定一组自由段射程角序列{β_i}；

其中，自由段射程角β为从主动段关机点起算的标准轨道射程角，对于所有(任何一个)真实轨道的自由段射程角，均需投影至标准轨道面来确定；自由段射程角间隔是根据计算精度与计算速度的要求来确定的，但一般不宜取的过小，因为过小会大大增加计算量而降低计算实时性，本实例中自由段射程角间隔取0.01度；

第二步，纵向摄动方程和侧向摄动方程的获得，具体是：

在轨道柱坐标系中描述弹道自由段的运动状态，轨道柱坐标系的具体定义为：如图1所示，坐标原点O与一个假想的沿着标准椭圆轨道a运动的弹道导弹的质心固连，它可以看作是导弹真实质心在标准轨道面内的投影，r轴在标准轨道面b内由地心O_E指向弹道导弹的质心O，即标准径向；β轴在标准轨道面内且垂直于r轴，沿着β增加的方向为正，即标准周向；z轴垂直于标准轨道面，且与r轴、β轴构成右手系，即标准侧向；

在该惯性系中，轨道柱坐标系三个坐标轴的单位矢量表示为表达式1)：

其中，β为自由段标准射程角；

它们的微分分别为：

以自由段起点为时间原点，任一时刻t的真实地心矢量由该时刻的标准轨道地心距r和标准侧向位移z表示如表达式3)：

则其速度表示为表达式4)：

其中为r的导数；

定义沿三个坐标轴方向的标准径向速度V_r、标准周向速度V_β和标准侧向速度V_z分别为表达式5)：

则其加速度表示为表达式6)：

其中为V_z的导数；

在考虑地球非球型摄动加速度影响的情况下，运用牛顿第二定律即可建立起受力状态与运动状态之间的关系如表达式7)：

其中，μ为地球引力常数；η为标准侧向偏差角，且δa_r、δa_β和δa_z分别为摄动加速度投影至轨道柱坐标系的三分量；

将标准侧向偏差角η和标准侧向速度V_z均看作小量，略去二阶及以上小量，综合表达式7)即得到以自由段标准射程角β为自变量的弹道导弹轨道柱坐标系简化运动模型，详情如表达式8)：

其中，上标“'”表示自由段标准射程角微分；

对于地球非球形摄动而言，相应的摄动加速度表示为非球形摄动引力位的梯度，再将其投影至标准轨道坐标系，得到分量形式的表达式9)：

基于简化的自由段运动方程，基于弹道摄动思想(详见图2)将此非线性方程线性化，并得出其状态转移矩阵解析解，具体是：

首先，在一条标准弹道附近，将关于摄动偏差量的弹道导弹自由段非线性运动方程进行线性化，得到等角/等地心距/等时摄动状态方程如表达式10)：

其中，δX为偏差状态量，V为摄动项；

其次，根据线性系统基本理论，推导该线性时变系统状态转移矩阵的解析解，即等角摄动/等地心距摄动/等时摄动的偏差半解析解如表达式11)：

其中，δX(β₀)为初始状态偏差；

这里给出等角摄动/等地心距摄动/等时摄动偏差的定义：

等角摄动偏差是指在自由段标准射程角相等的情况下，真实弹道与标准弹道的状态参数的偏差；等地心距摄动偏差是指真实弹道地心距与标准弹道地心距相等的情况下状态参数的偏差。等地心距摄动偏差更为适用于弹道摄动的计算，等地心距摄动偏差与等角摄动偏差的关系如图3所示；等时摄动偏差是指在自由段飞行时间相等的情况下，真实弹道与标准弹道的状态参数的偏差；

统一地记等角摄动/等地心距摄动/等时摄动的偏差状态Y为表达式12)：

Y＝[y₁ y₂ y₃ y₄ y₅ y₆]^Τ 12)，

以标准弹道为参考弹道，对简化的弹道导弹轨道柱坐标系自由段运动模型进行线性化，得等角摄动方程如表达式13)：

其中，

同样地，等地心距摄动状态方程如表达式14)：

其中：

等时摄动状态方程如表达式15)：

其中：

上述三种线性系统的系数矩阵并非常数，即为时变系统，通过对变量进行一定的转换获得便于求解状态转移矩阵解析表达式的形式，变化规则如下：

令

则，等角摄动方程通过适当的变量替换后，转变为表达式17)：

同时，给定初始条件

x_i(β₀)＝x_i(β₀)i＝1,2,…,6 18)，

显然，表达式17)中的第三式相对于其余五式是独立的，可以直接求解得到表达式19)：x₃(β)＝x₃(β₀) 19)，

将表达式17)中的第三式代入第一式中即可得到表达式20)：

将表达式20)代入表达式17)中的第二式可得表达式21)：

侧向摄动状态方程与纵向摄动状态方程相互解耦，因此可把表达式17)中的第六式代入第五式将侧向偏差参数悉数解出如表达式22)：

最后，将式表达式19)和表达式21)代入表达式17)中的第四式，完成纵向偏差参数的求解如表达式23)：

式中，

其中，f、r、p、h和e分别为真近点角、地心距、半通径、单位质量动量矩和偏心率；

经过上述求解流程和策略，实现了对新变量状态转移矩阵解析解的求解，将此状态转移矩阵整理为矩阵形式如表达式24)：

Φ(β,β₀)＝[Φ_ij(β,β₀)]i,j＝1,2,…,6 24)，

其中，

Φ₁₁(β,β₀)＝cos(β-β₀)；

Φ₁₂(β,β₀)＝-sin(β-β₀)；

Φ₂₁(β,β₀)＝sin(β-β₀)；

Φ₂₂(β,β₀)＝cos(β-β₀)；

Φ₃₃(β,β₀)＝1；

Φ₄₄(β,β₀)＝1；

Φ₅₅(β,β₀)＝cos(β-β₀)；

Φ₅₆(β,β₀)＝-sin(β-β₀)；

Φ₆₅(β,β₀)＝sin(β-β₀)；

Φ₆₆(β,β₀)＝cos(β-β₀)；

其余未列出各项均为0；

根据上述相同的解算思路，矩阵形式的等地心距摄动状态转移矩阵解析解可表示为表达式25)：

λ_r(β,β₀)＝[λ_rij(β,β₀)]i,j＝1,2,…,6 25)，

其中：

λ_r22(β,β₀)＝1；

λ_r44(β,β₀)＝1；

其余未列出各项均为0；

根据上述相同的解算思路，矩阵形式的等时摄动状态转移矩阵解析解表示为表达式26)：

λ_t(β,β₀)＝[λ_tij(β,β₀)]i,j＝1,2,…,6 26)，

其中：

λ_t44(β,β₀)＝1；

式中，

其余未列出各项均为0；

根据等角摄动/等地心距摄动/等时摄动模型，当自由段标准射程角为β时，等角摄动/等地心距摄动/等时摄动的过程偏差表示为表达式27)：

其中，U(ξ)为等角摄动/等地心距摄动/等时过程摄动项，具体如表达式28)：

将表达式28)代入表达式27)可得纵向摄动方程如表达式29)以及侧向摄动方程如表达式30)：

其中，μ为地球引力常数，p为标准椭圆弹道的半通径，r为地心距，y_i i＝1…4为纵向摄动偏差，δa_r和δa_β为摄动加速度投影在径向和周向的分量，λ为等角摄动/等地心距摄动/等时摄动状态转移矩阵解析解；

其中，μ为地球引力常数，p为标准椭圆弹道的半通径，r为地心距，y_i i＝5,6为侧向摄动偏差，δa_r为摄动加速度投影在径向的分量，λ为等角/等地心距/等时摄动状态转移矩阵解析解；

第三步，经过高阶项的最小二乘拟合后获得高阶摄动项的解析表达式，具体过程如下：

设观测方程为表达式31)：

z＝ha+ε 31)，

其中，z为观测向量，即以一定自由段标准射程角间隔积分解算半解析的侧向摄动方程所得到的一组侧向偏差角；h为系数矩阵，由线性拟合公式确定；a为待定参数向量；ε为随机误差向量；

设侧向偏差角描述为以自由段标准射程角β为自变量的简单函数的线性组合采用一个简单而可靠的经验拟合公式，如表达式32)：

该经验拟合公式可以看成对侧向偏差角的平方η²关于自由段标准射程角β的傅里叶级数形式，理论上级数的维数n达到无穷大时，以实现准确的描述，但在实际应用中往往截断到一定维数n，以满足运算精度的要求；所以，为了找到一个在计算量与计算精度之间达到平衡的拟合公式，需要确定一个合适的拟合阶数n，任意给定一组仿真条件下，对轨道倾角进行全方位循环，以比较各阶拟合公式的相对误差；从图4中以看出，5阶拟合公式显然在拟合精度上略低，而7阶与9阶拟合公式在精度上不相上下，相对误差均控制在0.2％以内，且7阶拟合公式的计算量明显优于9阶拟合公式，因此选择7阶拟合多项式作为拟合公式，即得线性拟合公式如表达式33)：

η²＝a₀+a₁sinβ+a₂cosβ+a₃sin2β+a₄cos2β+a₅sin3β+a₆cos3β 33)，

线性拟合公式确定之后，待定参数向量a与系数矩阵h也就随之确定，详见表达式34)和表达式37)，结合观测向量z与随机误差向量ε的表达式详见表达式35)和表达式36)，即可获得完整的观测方程，详情如下：

记a＝[a₀ a₁ a₂ a₃ a₄ a₅ a₆]^Τ 34)；

ε＝[ε₁ ε₂ … ε_n]^Τ 36)；

且

则根据最小二乘法，可以得到待定参数向量的线性无偏最优估计表达式38)：

a＝(h^Th)^-1h^Tz 38)；

从而获得高阶项的最小二乘估计式表达式39)：

此公式即为高阶摄动项的解析表达式；

第四步，对状态参数进行高阶偏差修正，具体如下：

当自由段标准射程角为β时，等地心距摄动高阶修正的过程偏差表示为表达式40)：

其中，W(β₀)为等地心距过程摄动项，△U_r(ξ)为高阶修正项，具体如表达式41)和表达式42)：

U_r(ξ)＝[△u_r 0 0 0 0 0]^Τ 42)，

高阶摄动项的最小二乘修正法的提出仅仅是为了修正径向参数的，对其它参数的修正几乎没有效果，因而忽略对其它状态偏差量的影响，仅考虑地心距偏差项的高阶修正即可，详见表达式43)：

将等地心距摄动状态转移矩阵解析解中λ_r21(β,ξ)的表达式代入表达式43)，即可得出高阶修正的等地心摄动纵向偏差参数如表达式44)：

其中：

同样可得出高阶修正的等角摄动纵向偏差参数详见表达式45)：

其中，u_i(β)如前等地心摄动纵向偏差参数表达式中的取值。

因等时摄动模型的高阶修正项的解析表达式不存在或者导不出，因此未列出。

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

Claims (10)

1.一种地球非球形摄动作用下自由段弹道误差传播的分析方法，其特征在于，包括以下步骤：

步骤一：选定自由段射程角序列，具体是：根据自由段射程角间隔Δβ和自由段射程角β，确定一组自由段射程角序列{β_i}；

步骤二：纵向摄动方程和侧向摄动方程的获得，具体是：建立以自由段标准射程角β为自变量的弹道导弹轨道柱坐标系简化运动模型，将简化的自由段运动方程线性化，并导出摄动模型和其状态转移矩阵解析解，包括：a、在一条标准弹道附近，将关于摄动偏差量的弹道导弹自由段非线性运动方程进行线性化，得到摄动模型下的摄动状态方程；b、根据线性系统基本理论，推导该线性时变系统状态转移矩阵的解析解；c、通过状态转移矩阵的解析解得到纵向摄动方程和侧向摄动方程；

步骤三：经过高阶项的最小二乘拟合后获得高阶摄动项的解析表达式，具体是：首先，设定观测方程和选定7阶拟合多项式作为拟合公式，其次，根据最小二乘法得到待定参数向量的线性无偏最优估计表达式；最后，获得高阶项的最小二乘估计式表达式，即为高阶摄动项的解析表达式；

步骤四：对状态参数进行高阶偏差修正，具体是：首先，取定当自由段标准射程角为β时高阶修正的过程偏差表达式；其次，对地心距偏差项进行高阶修正；最后，得到高阶修正的纵向偏差参数表达式。

2.根据权利要求1所述的地球非球形摄动作用下自由段弹道误差传播的分析方法，其特征在于，自由段射程角β为从主动段关机点起算的标准轨道射程角，所有真实轨道的自由段射程角均需投影至标准轨道面来进行确定；所述自由段射程角间隔Δβ为0.01度。

3.根据权利要求1所述的地球非球形摄动作用下自由段弹道误差传播的分析方法，其特征在于，所述步骤二中自由段运动方程的简化过程如下：

在轨道柱坐标系中，三个坐标轴的单位矢量表示为表达式1)：

其中，β为自由段标准射程角；表示坐标系O-XYZ三个坐标轴方向的单位矢量，表示坐标系O-rβz三个坐标轴方向的单位矢量；

三个坐标轴的单位矢量的微分分别为：

以自由段起点为时间原点，任一时刻t的真实地心矢量由该时刻的标准轨道地心距r和标准侧向位移z表示如表达式3)：

则其速度表示为表达式4)：

其中为r的导数；

定义沿三个坐标轴方向的标准径向速度V_r、标准周向速度V_β和标准侧向速度V_z分别为表达式5)：

则其加速度表示为表达式6)：

其中为V_z的导数；

在地球非球型摄动加速度影响的情况下，运用牛顿第二定律建立起受力状态与运动状态之间的关系如表达式7)：

其中，μ为地球引力常数；η为标准侧向偏差角，且δa_r、δa_β和δa_z分别为摄动加速度投影至轨道柱坐标系的三分量；

将标准侧向偏差角η和标准侧向速度V_z均看作小量，略去二阶及以上小量，得到以自由段标准射程角β为自变量的弹道导弹轨道柱坐标系简化运动模型，详情如表达式8)：

其中，上标“'”表示对自由段标准射程角取微分；

对于地球非球形摄动，相应的摄动加速度表示为非球形摄动引力位的梯度，再将其投影至标准轨道坐标系，得到分量形式的表达式9)：

其中，T表示地球非球型引力位。

4.根据权利要求3所述的地球非球形摄动作用下自由段弹道误差传播的分析方法，其特征在于，所述步骤二中的摄动状态方程具体为表达式10)：

其中，δX为偏差状态量，V为摄动项。

5.根据权利要求4所述的地球非球形摄动作用下自由段弹道误差传播的分析方法，其特征在于，所述步骤二中的摄动模型为等角摄动模型、等地心摄动模型以及等时摄动模型，三种摄动模型的偏差状态Y为表达式12)：

Y＝[y₁ y₂ y₃ y₄ y₅ y₆]^Τ 12)，

以标准弹道为参考弹道，对简化的弹道导弹轨道柱坐标系自由段运动模型进行线性化，得三种摄动模型下的摄动状态方程如下：

等角摄动状态方程如表达式13)：

其中，

等地心距摄动状态方程如表达式14)：

其中：

等时摄动状态方程如表达式15)：

其中：

6.根据权利要求5所述的地球非球形摄动作用下自由段弹道误差传播的分析方法，其特征在于，所述等角摄动状态转移矩阵的解析解的推导过程如下：

令

其中，u表示地心距的倒数；

则等角摄动状态方程转变为表达式17)：

同时，给定初始条件

x_i(β₀)＝x_i(β₀)i＝1,2,…,6 18)，

将表达式17)中的第三式求解得到表达式19)：

x₃(β)＝x₃(β₀) 19)，

将表达式17)中的第三式代入第一式中即得到表达式20)：

将表达式20)代入表达式17)中的第二式得表达式21)：

把表达式17)中的第六式代入第五式将侧向偏差参数解出如表达式22)：

将式表达式19)和表达式21)代入表达式17)中的第四式，完成纵向偏差参数的求解如表达式23)：

式中，

其中，f、r、p、h和e分别为真近点角、地心距、半通径、单位质量动量矩和偏心率；

将纵向偏差参数求解表达式23)整理为矩阵形式如表达式24)：

Φ(β,β₀)＝[Φ_ij(β,β₀)]i＝1,2,…,6；j＝1,2,…,6 24)，

其中，

Φ₁₁(β,β₀)＝cos(β-β₀)；

Φ₁₂(β,β₀)＝-sin(β-β₀)；

Φ₂₁(β,β₀)＝sin(β-β₀)；

Φ₂₂(β,β₀)＝cos(β-β₀)；

Φ₃₃(β,β₀)＝1；

Φ₄₄(β,β₀)＝1；

Φ₅₅(β,β₀)＝cos(β-β₀)；

Φ₅₆(β,β₀)＝-sin(β-β₀)；

Φ₆₅(β,β₀)＝sin(β-β₀)；

Φ₆₆(β,β₀)＝cos(β-β₀)；

其余未列出各项均为0。

7.根据权利要求6所述的地球非球形摄动作用下自由段弹道误差传播的分析方法，其特征在于，所述步骤二中：

矩阵形式的等地心距摄动状态转移矩阵解析解表示为表达式25)：

λ_r(β,β₀)＝[λ_rij(β,β₀)]其中：i＝1,2,…,6；j＝1,2,…,6 25)，

其中：

λ_r22(β,β₀)＝1；

λ_r44(β,β₀)＝1；

其余未列出各项均为0；

矩阵形式的等时摄动状态转移矩阵解析解表示为表达式26)：

λ_t(β,β₀)＝[λ_tij(β,β₀)]i＝1,2,…,6；j＝1,2,…,6 26)，

其中：

λ_t44(β,β₀)＝1；

式中，

其余未列出各项均为0。

8.根据权利要求7所述的地球非球形摄动作用下自由段弹道误差传播的分析方法，其特征在于，当自由段标准射程角为β时，等角摄动、等地心距摄动以及等时摄动三种摄动模型的过程偏差均表示为表达式27)：

其中，U(ξ)为等角摄动/等地心距摄动/等时过程摄动项，具体如表达式28)：

其中，r_ξ表示地心角为ξ时对应的地心距；

将表达式28)代入表达式27)得纵向摄动方程如表达式29)以及侧向摄动方程如表达式30)：

其中，μ为地球引力常数，p为标准椭圆弹道的半通径，r为地心距；y_i为纵向摄动偏差，其中i＝1…4；δa_r和δa_β为摄动加速度投影在径向和周向的分量，λ为等角摄动/等地心距摄动/等时摄动状态转移矩阵解析解；

其中，μ为地球引力常数，p为标准椭圆弹道的半通径，r为地心距；y_i为纵向摄动偏差，其中i＝5,6；δa_r为摄动加速度投影在径向的分量，λ为等角/等地心距/等时摄动状态转移矩阵解析解。

9.根据权利要求8所述的地球非球形摄动作用下自由段弹道误差传播的分析方法，其特征在于，所述步骤三中：

所述观测方程为表达式31)：

z＝ha+ε 31)，

其中，z为观测向量，即以一定自由段标准射程角间隔积分解算半解析的侧向摄动方程所得到的一组侧向偏差角；h为系数矩阵，由线性拟合公式确定；a为待定参数向量；ε为随机误差向量；

所述7阶拟合多项式作为拟合公式如表达式33)：

η²＝a₀+a₁sinβ+a₂cosβ+a₃sin2β+a₄cos2β+a₅sin3β+a₆cos3β 33)，

待定参数向量a与系数矩阵h详见表达式34)和表达式37)，观测向量z与随机误差向量ε的表达式详见表达式35)和表达式36)：

a＝[a₀ a₁ a₂ a₃ a₄ a₅ a₆]^Τ 34)；

ε＝[ε₁ ε₂ … ε_n]^Τ 36)；

且

所述待定参数向量的线性无偏最优估计表达式如表达式38)：

a＝(h^Th)^-1h^Tz 38)；

所述高阶项的最小二乘估计式表达式如表达式39)：

10.根据权利要求9所述的地球非球形摄动作用下自由段弹道误差传播的分析方法，其特征在于，所述步骤四中对等地心距摄动的状态参数进行高阶偏差修正具体过程如下：

当自由段标准射程角为β时，等地心距摄动高阶修正的过程偏差表示为表达式40)：

其中，W(ξ)为等地心距过程摄动项，U_r(ξ)为高阶修正项，具体如表达式41)和表达式42)：

U_r(ξ)＝[Δu_r(ξ) 0 0 0 0 0]^Τ 42)；

高阶摄动项的最小二乘修正法是为了修正径向参数，因此，地心距偏差项的高阶修正为表达式43)：

将等地心距摄动状态转移矩阵解析解中λ_r21(β,ξ)的表达式代入表达式43)即得高阶修正的等地心摄动纵向偏差参数如表达式44)：

其中：

u₀(β)＝1-cosβ；

所述步骤四中高阶修正的等角摄动纵向偏差参数详见表达式45)：

其中：