CN110044210B - 考虑任意阶地球非球型引力摄动的闭路制导在线补偿方法 - Google Patents

Abstract

本发明提供一种考虑任意阶地球非球型引力摄动的闭路制导在线补偿方法，包括以下步骤：第一步、沿弹道扰动引力重构模型离线构建；第二步、在弹上每一个制导周期内，基于二体动力学模型计算需要速度v_R；第三步、以弹上实时预测的关机点状态为输入条件，采用考虑地球非球型引力摄动的弹道偏差解析预报模型快速计算弹道终端状态偏差△X_c；第四步、根据弹道终端状态偏差△X_c解析计算需要速度修正量△v_R；第五步、求解考虑地球非球型引力摄动的控制量U。本发明基于弹道导弹自由段弹道偏差解析预报模型，提出了一种全新的闭路制导在线补偿方法，为进一步提升我国战略导弹制导方法的精度提供方法支撑。

Description

考虑任意阶地球非球型引力摄动的闭路制导在线补偿方法

技术领域

本发明涉及飞行动力学技术领域，具体涉及一种考虑任意阶地球非球型引力摄动的闭路制导在线补偿方法。

背景技术

弹上制导系统的性能是影响战略导弹作战效能的关键因素。不断提升导弹制导系统的精度、对弹道偏差的适应性是应对当前战略导弹快速机动发射及高精度命中需求的根本途径。闭路制导是根据目标数据和导弹的实时运动参数，按控制泛函(如射程偏差、需要速度或轨道运行周期等)的函数表达式进行实时计算的制导方法。该方法具有较高的适应性和鲁棒性，已成为当前战略导弹制导方法的首要选择。但闭路制导算法通常比较复杂，因此在进行制导指令的计算时很难顾及地球高阶非球型引力的影响，而基于弹道诸元修正的补偿思路对大范围弹道偏差的适应性较差，而且会增加战略导弹的发射准备时间，进而降低其射前生存能力。

针对当前存在的问题的不足，设计一种考虑任意阶地球非球型引力摄动的闭路制导在线补偿方法具有重要意义。

发明内容

本发明的目的在于提供一种考虑任意阶地球非球型引力摄动的闭路制导在线补偿方法，本发明基于弹道导弹自由段弹道偏差解析预报模型，提出了一种全新的闭路制导在线补偿方法，为进一步提升我国战略导弹制导方法的精度提供方法支撑。具体技术方案如下：

一种考虑任意阶地球非球型引力摄动的闭路制导在线补偿方法，包括以下步骤：

第一步、沿弹道扰动引力重构模型离线构建；第二步、在弹上每一个制导周期内，基于二体动力学模型计算需要速度v_R；第三步、以弹上实时预测的关机点状态为输入条件，采用考虑地球非球型引力摄动的弹道偏差解析预报模型快速计算弹道终端状态偏差△X_c；第四步、根据弹道终端状态偏差△X_c解析计算需要速度修正量△v_R；第五步、求解考虑地球非球型引力摄动的控制量U。

以上技术方案中优选的，沿弹道扰动引力重构模型离线构建具体是：根据导弹飞行任务设计一条标准弹道，以标准弹道关机点参数为基准，在导弹发射前构建扰动引力重构模型，并将每个网格节点的位置坐标矢量及扰动引力矢量信息存储在弹上。

以上技术方案中优选的，基于弹上导航系统解算出的导弹当前位置矢量r₀、目标点位置矢量r_F、导弹总飞行时间t_total及当前时刻t，采用Lambert算法计算当前点处的需要速度v_R，其中导弹总飞行时间为提前设计的固定值，Lambert问题求解时输入的时间项为导弹总飞行时间减去当前时刻，即t_total-t；Lambert问题求解时，需要速度

和飞行时间

采用表达式1)和表达式2)进行计算：

以上技术方案中优选的，地球非球型引力引起的弹道终端偏差△X_c采用表达式3)进行计算：

以上技术方案中优选的，地球J₂项引力引起的弹道偏差△X_J2的计算方法包括J₂项引力矢量的分解以及自由段弹道偏差解析预报模型推导；

J₂项引力矢量的分解得到不同坐标轴方向上的摄动力如表达式3-1-8)：

自由段弹道偏差解析预报模型推导具体是：

根据状态空间摄动理论，得到考虑J₂项影响的自由段弹道偏差解析解如表达式3-1-9)-3-1-14)：

则有表达式3-1-15)：

△X_J2＝[△v_r(f),△r(f),△v_β(f),△t(f),△v_z(f),△z(f)]^T 3-1-15)。

以上技术方案中优选的，J₂项引力矢量的分解，具体是：

令U₂表示地球J₂项引力势，在地球惯性系中J₂项引力势为表达式3-1-1)：

根据球面三角学将表达式1)的自变量由地心纬度转换为真近点角f；取N为天球坐标系的北天极，大圆弧

为由导弹关机点参数确定的标准二体弹道在球面上的投影，曲线AB表示导弹受摄运动轨迹在球面上的投影，直线OP垂直于平面OAB^*，α_A、△f分别为平面POC与POA的二面角和平面POA与POB^*的二面角，

分别为A点的纬度和P点的纬度，λ_A和λ_P分别为A点和P点的经度，则有：

在球面三角形ANP中，得表达式3-1-2)至3-1-4)：

在球面三角形BPN中，得表达式3-1-5)：

将表达式3-1-5)带入表达式3-1-1)，即得到J₂项引力势关于真近点角的函数为表达式3-1-6)：

分别求U₂(f)关于r、f、σ的偏导数，即得到J₂项引力矢量在轨道柱坐标系中的表达式为表达式3-1-7)：

基于二体标准弹道进行计算，此时σ＝0，表达式3-1-7)中的如下系数退化为零，即

将不同坐标轴方向上的摄动力统一用表达式3-1-8)表示：

以上技术方案中优选的，由扰动引力引起的弹道偏差△X_δg的计算方法包括高阶扰动引力矢量分解和状态偏差解析预报模型推导；

任意阶扰动引力矢量的分解包括以下过程：

先得到P点处扰动引力矢量δg关于真近点角的函数表达式如表达式3-2-4)：

再获得扰动引力三分量的表达式为表达式3-2-5)：

最后得到表达式3-2-6)：

状态偏差解析预报模型推导具体包括以下步骤：

根据状态空间摄动理论，获得考虑高阶扰动引力影响的自由段弹道偏差解析解如表达式3-2-7)至3-2-12)：

则有表达式3-2-13)：

△X_δg＝[△v′_r(f),△r′(f),△v′_β(f),△t′(f),△v′_z(f),△z′(f)]^T 3-2-13)。

以上技术方案中优选的，P点处扰动引力矢量δg关于真近点角的函数表达式3-2-4)获得的具体过程如下：

根据沿飞行弹道的扰动引力重构模型，弹道上任意一点的扰动引力矢量表示为八面体网格八个节点扰动引力矢量的加权和，即为表达式3-2-1)：

其中：插值核函数

满足

式中，ξ_k、

和η_k为第k个节点的位置坐标，而ξ′_k、

和η′_k分别为与该节点相邻的三个节点对应坐标轴方向上的坐标值，δg为弹道上任意一点的扰动引力矢量，δg_k为有限元网格八个节点的扰动引力矢量：

弹道上任意点P与轨道坐标系

之间的几何关系图中，诸元二体弹道是基于标准关机点K_f确定的二体弹道，用于在导弹发射前进行扰动引力重构模型的构建；制导二体弹道是由导弹实际关机点K确定的二体弹道；实际弹道为导弹在地球非球型引力等摄动力作用下的真实飞行弹道；ψ₀和φ₀为地心角；P′为制导二体弹道上与P点对应的点，设P′在

系中的位置矢量为

则有表达式3-2-2)：

式中：

为坐标系

原点对应的地心距；Μ_z(·)和Μ_y(·)分别表示绕z轴和y轴变换的方向余弦矩阵，且有

表示真近点角为f时对应的标准二体弹道地心距，即

△β、ψ₀和φ₀均为矢量夹角，△f_i表示当前有限元网格截得的诸元标准二体弹道地心角；ψ₀和φ₀由K_f和K之间的几何关系精确计算出来，具体计算过程如下：

设K_f点处的位置矢量和速度矢量分别为r_k和v_k，则诸元二体轨道平面的动量矩矢量h_k为h_k＝r_k×v_k；再设K点处的位置矢量为r₀，则有

令矢量r_k和r₀之间的夹角为θ₀，可得

根据ψ₀、θ₀和φ₀之间的几何关系，可得

△β满足△β＝f-f_K，其中f和f_K分别为点P′和K在制导二体弹道上对应的真近点角；

将表达式3-2-2)展开有表达式3-2-3)：

式中：r为地球半径，p₁-p₆如下：

将表达式3-2-3)带入表达式3-2-1)得到P点处扰动引力矢量δg关于真近点角的函数表达式如表达式3-2-4)：

其中，α_i均为常矢量系数，i＝0,1,…,9。

以上技术方案中优选的，增益速度修正量δ△v采用表达式7)进行计算：

δ△v＝[δ△v_r0,δ△v_β0,δ△v_z0]^T 7)；

其中：δ△v_r0、δ△v_β0和δ△v_z0分别表示沿二体轨道径向、切向和侧向的初始速度增量，则有

式中：

以上技术方案中优选的，考虑地球非球型引力摄动的控制量U具体是根据修正后的增益速度△v′_R＝△v_R+δ△v计算得到，具体是：

根据发射惯性系与轨道柱坐标系的关系得到表达式8)：

式中：δ△v′为在发射惯性系中表示的增益速度修正量，

为轨道柱坐标系到地心惯性系的方向余弦阵，

地心惯性系到发射惯性系的方向余弦阵，具体是：

设修正后的增益速度矢量为△v′_R＝[△v_x △v_y △v_z]^T，则得控制量

为表达式9)：

式中：

和ψ_c分别为当前导弹实际俯仰角和偏航角。

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

附图说明

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

图1是实施例中考虑地球非球型引摄动在线补偿的闭路制导系统的结构示意图；

图2是实施例Lambert问题的求解流程示意图；

图3是实施例中地球J₂项引力引起的弹道偏差获得时的球面角关系图；

图4是实施例中地球J₂项引力引起的弹道偏差获得时的轨道柱坐标系示意图；

图5是实施例中由扰动引力引起的弹道偏差获得时弹道任意点位置坐标转换图；

图6是实施例中俯仰角随制导过程的变化曲线图；

图7是实施例中偏航角随制导过程的变化曲线图；

图8是实施例中X方向增益速度随制导过程的变化曲线图；

图9是实施例中Y方向增益速度随制导过程的变化曲线图；

图10是实施例中标称条件下补偿J₂项和扰动引力条件下的射程偏差补偿效果图；

图11是实施例中标称条件下补偿J₂项和扰动引力条件下的横向偏差补偿效果图；

图12是实施例中标称条件下仅补偿J₂项引力条件下的射程偏差补偿效果图；

图13是实施例中标称条件下仅补偿J₂项引力条件下的横向偏差补偿效果图；

图14是实施例中大偏差条件下两种方法射程补偿残差图；

图15是实施例中大偏差条件下两种方法横向补偿残差图。

具体实施方式

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

实施例：

一种考虑任意阶地球非球型引力摄动的闭路制导在线补偿方法，详见图1，包括以下步骤：

第一步、沿弹道扰动引力重构模型离线构建，具体是：根据导弹飞行任务设计一条标准弹道，以标准弹道关机点参数为基准，在导弹发射前构建扰动引力重构模型(可参照专利“沿飞行弹道的扰动引力重构模型”中描述的方法)，并将每个网格节点的位置坐标矢量及扰动引力矢量信息存储在弹上。

第二步、在弹上每一个制导周期内，基于二体动力学模型计算需要速度v_R，具体是：基于弹上导航系统解算出的导弹当前位置矢量r、目标点位置矢量r_c、导弹总飞行时间t_total及当前时刻t，采用Lambert算法计算当前点处的需要速度v_R，其中导弹总飞行时间为提前设计的固定值，Lambert问题求解时输入的时间项为导弹总飞行时间减去当前时刻，即t_total-t。

图2所示为导弹飞行时间固定条件下的Lambert问题求解流程。针对弹道导弹的运动特性，速度倾角γ存在上界和下界，且其值满足如下表达式：

通过在区间(γ_min,γ_max)内迭代速度倾角，便可实现Lambert问题的快速求解。迭代初值γ₀一般取为γ₀＝(γ_min+γ_max)/2。每一次迭代中，需要速度

和飞行时间

采用表达式1)和表达式2)进行计算：

第三步、以弹上实时预测的关机点状态为输入条件，采用考虑地球非球型引力摄动的弹道偏差解析预报模型快速计算弹道终端状态偏差△X_c，具体是：

地球非球型引力引起的弹道终端偏差△X_c采用表达式3)进行计算：

地球J₂项引力引起的弹道偏差△X_J2的计算方法包括J₂项引力矢量的分解以及自由段弹道偏差解析预报模型推导；

J₂项引力矢量的分解具体如下：

令U₂表示地球J₂项引力势，在地球惯性系中J₂项引力势为表达式3-1-1)：

根据球面三角学将表达式1)的自变量由地心纬度转换为真近点角f；如图3，取N为天球坐标系的北天极，大圆弧

为由导弹关机点参数确定的标准二体弹道在球面上的投影，曲线AB表示导弹受摄运动轨迹在球面上的投影，直线OP垂直于平面OAB^*，α_A、△f分别为平面POC与POA的二面角和平面POA与POB^*的二面角，

分别为A点的纬度和P点的纬度，λ_A和λ_P分别为A点和P点的经度，则有：

在球面三角形ANP中，得表达式3-1-2)至3-1-4)：

在球面三角形BPN中，得表达式3-1-5)：

将表达式3-1-5)带入表达式3-1-1)，即得到J₂项引力势关于真近点角的函数为表达式3-1-6)：

分别求U₂(f)关于r、f、σ的偏导数，即得到J₂项引力矢量在轨道柱坐标系(轨道柱坐标系的定义如图4所示)中的表达式为表达式3-1-7)：

基于状态空间摄动法的弹道误差传播解析解推导过程中，摄动力并非基于当前导弹真实位置进行计算，而是基于二体标准弹道进行计算，此时σ＝0，表达式3-1-7)中的如下系数退化为零，即

将不同坐标轴方向上的摄动力统一用表达式3-1-8)表示：

自由段弹道偏差解析预报模型推导具体是：

根据状态空间摄动理论，导弹自由飞行段弹道偏差的积分求解表示式如下：

式中：h为二体弹道平面对应的动量矩矢量的模；

表示真近点角为ξ时对应的标准二体弹道地心距，即

p表示二体弹道的半通径，e表示二体弹道的偏心率；

μ为地球引力常数；

将表达式3-1-8)代入导弹自由飞行段弹道偏差的积分求解表示式，并积分即得每一项偏差的完整解析表达式3-1-9)至表达式3-1-14)：

则有表达式3-1-15)：

△X_J2＝[△v_r(f),△r(f),△v_β(f),△t(f),△v_z(f),△z(f)]^T 3-1-15)。

由扰动引力引起的弹道偏差△X_δg的计算方法包括高阶扰动引力矢量分解和状态偏差解析预报模型推导；

任意阶扰动引力矢量的分解包括以下过程：

先得到P点处扰动引力矢量δg关于真近点角的函数表达式，具体过程如下：

根据沿飞行弹道的扰动引力重构模型，弹道上任意一点的扰动引力矢量表示为八面体网格八个节点扰动引力矢量的加权和，即为表达式3-2-1)：

其中：插值核函数

满足

式中，ξ_k、

和η_k为第k个节点的位置坐标，而ξ′_k、

和η′_k分别为与该节点相邻的三个节点对应坐标轴方向上的坐标值，δg为弹道上任意一点的扰动引力矢量，δg_k为有限元网格八个节点的扰动引力矢量：

图5所示为弹道上任意点P与轨道坐标系

之间的几何关系图中，其中诸元二体弹道是基于标准关机点K_f确定的二体弹道，用于在导弹发射前进行扰动引力重构模型的构建；制导二体弹道是由导弹实际关机点K确定的二体弹道；实际弹道为导弹在地球非球型引力等摄动力作用下的真实飞行弹道；ψ₀和φ₀为地心角；P′为制导二体弹道上与P点对应的点，设P′在

系中的位置矢量为

则根据图5所示的几何关系有表达式3-2-2)：

表达式中ψ₀和φ₀由K_f和K之间的几何关系精确计算出来，具体计算过程如下：

设K_f点处的位置矢量和速度矢量分别为r_k和v_k，则诸元二体轨道平面的动量矩矢量h_k为h_k＝r_k×v_k；再设K点处的位置矢量为r₀，则有

令矢量r_k和r₀之间的夹角为θ₀，可得

根据ψ₀、θ₀和φ₀之间的几何关系，可得

△β满足△β＝f-f_K，其中f和f_K分别为点P′和K在制导二体弹道上对应的真近点角；

将表达式3-2-2)展开有表达式3-2-3)：

将表达式3-2-3)带入表达式3-2-1)得到P点处扰动引力矢量δg关于真近点角的函数表达式如表达式3-2-4)：

扰动引力影响下自由段弹道偏差解析解的推导需要扰动引力三分量的表达式，令u^φ(φ＝r,β,z)表示标准二体弹道上任意点处扰动引力的三分量，

表示每个八面体网格节点上扰动引力的三分量，则有表达式3-2-5)：

最后得到表达式3-2-6)：

状态偏差解析预报模型推导具体包括以下步骤：

根据状态空间摄动理论，获得导弹自由飞行段弹道偏差的积分求解表示式如下：

将λ_j,s(f,ξ),j＝1,2,…,6,s＝1,3,5以及表达式3-2-6)代入导弹自由飞行段弹道偏差的积分求解表示式中积分得到考虑高阶扰动引力影响的自由段弹道偏差解析解如表达式3-2-7)至3-2-12)：

则有表达式3-2-13)：

△X_δg＝[△v′_r(f),△r′(f),△v′_β(f),△t′(f),△v′_z(f),△z′(f)]^T 3-2-13)。

第四步、根据弹道终端状态偏差△X_c解析计算需要速度修正量△v_R，具体是：修正增益速度的目的是消除由地球非球型引力摄动引起的弹道终端位置偏差，但由于地球存在自转，因此时间项偏差也会导致落点偏差。

假设

则根据状态空间摄动法，受摄情况下导弹状态偏差的计算公式为表达式4)：

式中：Φ(f,f₀)为线性微分方程的状态转移矩阵，f和f₀分别为初始时刻和终端时刻对应的二体轨道真近点角；U(ξ)表示摄动力，此处表示地球非球型引力。

显然，由第二步计算出的导弹终端状态偏差△X_c满足表达式5)：

由表达式5)可知，要使△X_f≡0，需要Φ(f,f₀)△X_f0＝-△X_c，表明可以通过施加一个初始状态偏差量来抵消由地球非球型引力引起的终端状态偏差。设沿二体轨道径向、切向和侧向的初始速度增量分别为δ△v_r0、δ△v_β0和δ△v_z0，则根据表达式6)可知：

式中，△v′_rc、△v′_βc、△v′_zc分别为由δ△v引起的弹道终端速度偏差，λ_i,j,i,j＝1,2,…,6为状态转移矩阵的成员变量。

求解表达式6)即可得增益速度修正量δ△v采用表达式7)进行计算：

δ△v＝[δ△v_r0,δ△v_β0,δ△v_z0]^T 7)；

其中：δ△v_r0、δ△v_β0和δ△v_z0分别表示沿二体轨道径向、切向和侧向的初始速度增量，则有

第五步、求解考虑地球非球型引力摄动的控制量U，具体是根据修正后的增益速度△v′_R＝△v_R+δ△v计算得到，详情如下：

易知δ△v是在轨道柱坐标系中描述的，控制量计算时需要将其转换至发射惯性系。根据发射惯性系与轨道柱坐标系的关系得到表达式8)：

式中：δ△v′为在发射惯性系中表示的增益速度修正量，

为轨道柱坐标系到地心惯性系的方向余弦阵，

地心惯性系到发射惯性系的方向余弦阵，具体是：

设修正后的增益速度矢量为△v′_R＝[△v_x △v_y △v_z]^T，则得控制量

为表达式9)：

式中

和ψ_c分别为当前导弹实际俯仰角和偏航角。

采用本实施例方法进行仿真中假定的战略导弹发射点大地经度为λ₀＝90°，大地纬度为

大地高度为100m，发射方位角由-180°遍历至180°。弹道导弹本体质量及发动机参数采用美国民兵-III洲际弹道导弹的模型参数，如表2所示。射程设置为12000km，飞行时间固定为2900秒。假设弹道导弹在第三级开始进行闭路制导。为了分析J₂项引力和扰动引力的影响量级，根据所考虑的自由段非球型引力模型将仿真分为三组：第一组仅考虑J₂项引力，第二组考虑72×72阶非球型引力，第三组考虑72×72阶扰动引力。

表2民兵3导弹发动机参数

级数	点火前质量(吨)	弹体直径(m)	发动机推力(kN)	发动机工作时间(s)
					第一级	34.5	1.67	912	61.6
第二级	11.82	1.32	270	65.2
					第三级	4.77	1.32	155	59.6

弹道导弹导弹从第三级开始闭路制导，开始时间为128秒。图6和图7所示分别为发射方位角为45度时，弹道导弹闭路制导过程中俯仰角和偏航角的变化曲线，其中实线为不考虑地球非球型引力影响补偿时的俯仰角和偏航角变化情况，虚线为考虑地球非球型引力影响在线补偿时俯仰角和偏航角变化情况。

图8和图9分别为发射方位角为45度时，弹道导弹闭路制导过程中x方向增益速度和y方向增益速度的变化曲线，其中实线为不考虑地球非球型引力影响补偿时的增益速度变化情况，虚线为考虑地球非球型引力影响在线补偿时的增益速度变化情况。

图10-图13分别为本发明所述方法与闭路制导地面诸元补偿方法在不同条件下的射程补偿残差和横向补偿残差随发射方位角的变化曲线，表3为两种方法在不同条件下的补偿残差统计结果，统计样本均为残差的绝对值。

表3标称条件下不同方法的射程补偿残差统计分析结果

图14-图15分别为本发明所述方法与闭路制导地面诸元补偿方法在考虑发动机推力偏差条件下的射程补偿残差和横向补偿残差随发射方位角的变化曲线，表4为两种方法在不同条件下的补偿残差统计结果，统计样本均为残差的绝对值。

表4大偏差条件下不同补偿方法计算精度统计分析结果

从仿真结果可以看出：

①在同时补偿J₂项引力和扰动引力的条件下，在线补偿方法的射程补偿残差最大值为25.0252米，补偿残差均值为10.9294米；横向补偿残差最大值为39.2784米，补偿残差均值为14.6647米。

②在仅补偿J₂项引力的条件下，在线补偿方法的射程补偿残差最大值为376.429米，补偿残差均值为342.623米；横向补偿残差最大值为291.382米，补偿残差均值为131.975米。显然，在闭路制导中忽略自由段扰动引力对需要速度的影响是导致该补偿残差的最主要原因；

③标称条件下，闭路制导地面诸元补偿精度比在线补偿方法稍差，但补偿残差也均在50米以内，精度较高。

④在大偏差条件下，地面诸元补偿方法的射程补偿残差最大值达到257.8639米，均值为224.8207米；横向补偿残差最大值为94.6340米，均值为36.0602米，精度较无偏差条件的情况低了一个量级。这主要是由于弹道变形会使得地球非球型引力对战略导弹自由飞段弹道的影响与标称条件下的影响产生偏差，因此地面计算的虚拟目标点并不能精确补偿实际地球非球型引力的影响，这使得地面诸元补偿方法对助推段弹道变形的适应能力较弱。

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

Claims (10)

1.一种考虑任意阶地球非球型引力摄动的闭路制导在线补偿方法，包括沿弹道扰动引力重构模型离线构建，其特征在于：还包括以下步骤：

在弹上每一个制导周期内，基于二体动力学模型计算需要速度v_R；

以弹上实时预测的关机点状态为输入条件，采用考虑地球非球型引力摄动的弹道偏差解析预报模型快速计算弹道终端状态偏差ΔX_c；

根据弹道终端状态偏差ΔX_c解析计算需要速度修正量Δv_R；

求解考虑地球非球型引力摄动的控制量U。

2.根据权利要求1所述的考虑任意阶地球非球型引力摄动的闭路制导在线补偿方法，其特征在于：沿弹道扰动引力重构模型离线构建具体是：根据导弹飞行任务设计一条标准弹道，以标准弹道关机点参数为基准，在导弹发射前构建扰动引力重构模型，并将每个网格节点的位置坐标矢量及扰动引力矢量信息存储在弹上。

3.根据权利要求1所述的考虑任意阶地球非球型引力摄动的闭路制导在线补偿方法，其特征在于：基于弹上导航系统解算出的导弹当前位置矢量r₀、目标点位置矢量r_F、导弹总飞行时间t_total及当前时刻t，采用Lambert算法计算当前点处的需要速度v_R，其中导弹总飞行时间为提前设计的固定值，Lambert问题求解时输入的时间项为导弹总飞行时间减去当前时刻，即t_total-t；Lambert问题求解时，需要速度

和飞行时间

采用表达式1)和表达式2)进行计算：

其中：速度倾角γ存在上界和下界，且其值满足如下表达式：

φ采用表达式

进行计算，r₀为导弹当前位置对应的地心距，μ为地球引力常数，r_F为导弹目标位置对应的地心距，λ取值为

γ_n采用公式

4.根据权利要求1所述的考虑任意阶地球非球型引力摄动的闭路制导在线补偿方法，其特征在于：地球非球型引力引起的弹道终端偏差ΔX_c采用表达式3)进行计算：

其中：ΔX_J2为地球J₂项引力引起的弹道偏差，ΔX_δg为高阶扰动引力引起的弹道偏差，Δv_r(f)、Δv_β(f)和Δv_z(f)分别为考虑J₂项影响的弹道状态偏差速度矢量在轨道柱坐标系中沿r轴、β轴和z轴方向的分量；Δr(f)和Δz(f)分别为考虑J₂项影响的弹道状态偏差位置矢量在轨道柱坐标系中沿r轴和z轴方向的分量；Δt(f)为考虑J₂项影响的实际飞行时间与标准二体弹道飞行时间之差，Δv′_r(f)、Δv′_β(f)和Δv′_z(f)分别为考虑高阶扰动引力影响的弹道状态偏差速度矢量在轨道柱坐标系中沿r轴、β轴和z轴方向的分量；Δr′(f)和Δz′(f)分别为考虑高阶扰动引力影响的弹道状态偏差位置矢量在轨道柱坐标系中沿r轴和z轴方向的分量；Δt′(f)为考虑高阶扰动引力影响的实际飞行时间与标准二体弹道飞行时间之差。

5.根据权利要求4所述的考虑任意阶地球非球型引力摄动的闭路制导在线补偿方法，其特征在于：地球J₂项引力引起的弹道偏差ΔX_J2的计算方法包括J₂项引力矢量的分解以及自由段弹道偏差解析预报模型推导；

J₂项引力矢量的分解得到不同坐标轴方向上的摄动力如表达式3-1-8)：

式中：δa_r、δa_β和δa_z分别表示J₂项引力矢量在轨道柱坐标系中r轴、β轴和z轴方向上的分量；r为地球半径；

s_r＝-3K,s_β＝K,s_z＝K；

μ为地球引力常数，a_e为地球赤道平均半径；p_i为常系数，i取0,1,2,3,4，具体如下：

为P点的纬度，σ为侧向角，α_A为A点在极点坐标系中的经度，f₀表示A点处的真近点角，即初始真近点角；

q₁和q₂如下：

自由段弹道偏差解析预报模型推导具体是：

根据状态空间摄动理论，得到考虑J₂项影响的自由段弹道偏差解析解如表达式3-1-9)-3-1-14)：

式中：

Λ_i,j如表1所示：

表1函数Λ_i,j的表达式统计表

则有表达式3-1-15)：

ΔX_J2＝[Δv_r(f),Δr(f),Δv_β(f),Δt(f),Δv_z(f),Δz(f)]^T 3-1-15)。

6.根据权利要求5所述的考虑任意阶地球非球型引力摄动的闭路制导在线补偿方法，其特征在于：J₂项引力矢量的分解，具体是：

令U₂表示地球J₂项引力势，在地球惯性系中J₂项引力势为表达式3-1-1)：

其中：μ为地球引力常数，a_e为地球赤道平均半径，J为常数且

r为地球半径，

为地球纬度；

根据球面三角学将表达式1)的自变量由地心纬度转换为真近点角f；取N为天球坐标系的北天极，大圆弧AB^*为由导弹关机点参数确定的标准二体弹道在球面上的投影，曲线AB表示导弹受摄运动轨迹在球面上的投影，直线OP垂直于平面OAB^*，α_A、Δf分别为平面POC与POA的二面角和平面POA与POB^*的二面角，

分别为A点的纬度和P点的纬度，λ_A和λ_P分别为A点和P点的经度，则有：

在球面三角形ANP中，得表达式3-1-2)至3-1-4)：

其中：γ为关机点A对应的方位角，α_A为A点在新极点坐标系中的经度；

在球面三角形BPN中，得表达式3-1-5)：

式中：σ为侧向角；f₀表示A点处的真近点角，即初始真近点角；f表示弹上任意时刻对应的真近点角，这里P即表示任意点；

将表达式3-1-5)带入表达式3-1-1)，即得到J₂项引力势关于真近点角的函数为表达式3-1-6)：

分别求U₂(f)关于r、f、σ的偏导数，即得到J₂项引力矢量在轨道柱坐标系中的表达式为表达式3-1-7)：

基于二体标准弹道进行计算，此时σ＝0，表达式3-1-7)中的如下系数退化为零，即

将不同坐标轴方向上的摄动力统一用表达式3-1-8)表示：

7.根据权利要求4所述的考虑任意阶地球非球型引力摄动的闭路制导在线补偿方法，其特征在于：由扰动引力引起的弹道偏差ΔX_δg的计算方法包括高阶扰动引力矢量分解和状态偏差解析预报模型推导；

任意阶扰动引力矢量的分解包括以下过程：

先得到P点处扰动引力矢量δg关于真近点角的函数表达式如表达式3-2-4)：

其中，α_i均为常矢量系数，i＝0,1,…,9；

表示真近点角为f时对应的标准二体弹道地心距，即

且p表示二体弹道的半通径，e表示二体轨道偏心率；f为点P′在制导二体弹道上对应的真近点角，P′为制导二体弹道上与P点对应的点；

再获得扰动引力三分量的表达式为表达式3-2-5)：

其中：u^φ(φ＝r,β,z)表示标准二体弹道上任意点处扰动引力的三分量；取变量n₁-n₈和α_i,k(i和k分别取0,1,2,3,4,5,6,7,8,9)如下：

α_6,i＝p₁p₄p₆+p₂p₃p₆+p₂p₄p₅；α_7,i＝p₁p₃p₆+p₁p₄p₅+p₂p₃p₅；

α_8,i＝p₂p₄p₆；α_9,i＝p₁p₃p₅；

ξ₁-ξ₈、

以及η₁-η₈为八面体网格八个节点在局部坐标系中的坐标位置；

最后得到表达式3-2-6)：

状态偏差解析预报模型推导具体包括以下步骤：

根据状态空间摄动理论，获得考虑高阶扰动引力影响的自由段弹道偏差解析解如表达式3-2-7)至3-2-12)：

式中：ε₁-ε₄为

P₁₀-P₁₉和Q₁₀-Q₁₉为

P₂₀-P₂₉和Q₂₀-Q₂₉为

P₃₀-P₃₉和Q₃₀-Q₃₉为

P₄₀-P₄₉和Q₄₀-Q₄₉为

P₅₀-P₅₉和Q₅₀-Q₅₉为

P₆₀-P₆₉和Q₆₀-Q₆₉为

所涉及的l^n,p,q(E)函数和κ^n,p,q(E)函数的解析表示式如下：

κ^2,0,1(E)＝a²(1-e²)(sinE-EcosE)；

a表示二体轨道半长轴，E表示偏近点角；

则有表达式3-2-13)：

ΔX_δg＝[Δv′_r(f),Δr′(f),Δv′_β(f),Δt′(f),Δv′_z(f),Δz′(f)]^T 3-2-13)。

8.根据权利要求7所述的考虑任意阶地球非球型引力摄动的闭路制导在线补偿方法，其特征在于：P点处扰动引力矢量δg关于真近点角的函数表达式3-2-4)获得的具体过程如下：

根据沿飞行弹道的扰动引力重构模型，弹道上任意一点的扰动引力矢量表示为八面体网格八个节点扰动引力矢量的加权和，即为表达式3-2-1)：

其中：插值核函数

满足

式中，ξ_k、

和η_k为第k个节点的位置坐标，而ξ′_k、

和η′_k分别为与该节点相邻的三个节点对应坐标轴方向上的坐标值，δg为弹道上任意一点的扰动引力矢量，δg_k为有限元网格八个节点的扰动引力矢量：

弹道上任意点P与轨道坐标系

之间的几何关系图中，诸元二体弹道是基于标准关机点K_f确定的二体弹道，用于在导弹发射前进行扰动引力重构模型的构建；制导二体弹道是由导弹实际关机点K确定的二体弹道；实际弹道为导弹在地球非球型引力等摄动力作用下的真实飞行弹道；ψ₀和φ₀为地心角；P′为制导二体弹道上与P点对应的点，设P′在

系中的位置矢量为

则有表达式3-2-2)：

式中：

为坐标系

原点对应的地心距；Μ_z(·)和Μ_y(·)分别表示绕z轴和y轴变换的方向余弦矩阵，且有

表示真近点角为f时对应的标准二体弹道地心距，即

Δβ、ψ₀和φ₀均为矢量夹角，Δf_i表示当前有限元网格截得的诸元标准二体弹道地心角；ψ₀和φ₀由K_f和K之间的几何关系精确计算出来，具体计算过程如下：

设K_f点处的位置矢量和速度矢量分别为r_k和v_k，则诸元二体轨道平面的动量矩矢量h_k为h_k＝r_k×v_k；再设K点处的位置矢量为r₀，则有

令矢量r_k和r₀之间的夹角为θ₀，可得

根据ψ₀、θ₀和φ₀之间的几何关系，可得

Δβ满足Δβ＝f-f_K，其中f和f_K分别为点P′和K在制导二体弹道上对应的真近点角；

将表达式3-2-2)展开有表达式3-2-3)：

式中：r为地球半径，p₁-p₆如下：

将表达式3-2-3)带入表达式3-2-1)得到P点处扰动引力矢量δg关于真近点角的函数表达式如表达式3-2-4)：

其中，α_i均为常矢量系数，i＝0,1,…,9。

9.根据权利要求5所述的考虑任意阶地球非球型引力摄动的闭路制导在线补偿方法，其特征在于：增益速度修正量δΔv采用表达式7)进行计算：

δΔv＝[δΔv_r0,δΔv_β0,δΔv_z0]^T 7)；

其中：δΔv_r0、δΔv_β0和δΔv_z0分别表示沿二体轨道径向、切向和侧向的初始速度增量，则有

式中：

10.根据权利要求9所述的考虑任意阶地球非球型引力摄动的闭路制导在线补偿方法，其特征在于：考虑地球非球型引力摄动的控制量U具体是根据修正后的增益速度Δv′_R＝Δv_R+δΔv计算得到，具体是：

根据发射惯性系与轨道柱坐标系的关系得到表达式8)：

式中：δΔv′为在发射惯性系中表示的增益速度修正量，

为轨道柱坐标系到地心惯性系的方向余弦阵，

地心惯性系到发射惯性系的方向余弦阵，具体是：

设修正后的增益速度矢量为Δv′_R＝[Δv_x Δv_y Δv_z]^T，则得控制量

为表达式9)：

式中：

和ψ_c分别为当前导弹实际俯仰角和偏航角；Δv_x、Δv_y和Δv_z为在发射惯性系中增益速度矢量在x轴、y轴和z轴方向上的分量。

Images