CN110046439A - 考虑高阶扰动引力影响的弹道偏差解析预报算法 - Google Patents

Abstract

本发明提供一种考虑高阶扰动引力影响的弹道偏差解析预报算法，包括高阶扰动引力矢量分解和状态偏差解析预报模型推导，任意阶扰动引力矢量的分解包括以下过程包括：先得到二体弹道上任意点P处扰动引力矢量δg关于真近点角的函数表达式，再获得扰动引力三分量的表达式，最后得到的表达式；状态偏差解析预报模型推导具体包括以下步骤：根据状态空间摄动理论，获得导弹自由飞行段弹道偏差的积分求解表示式，将代入并积分即得考虑扰动引力影响的自由段弹道偏差解析解的表达式。本发明方法的解算效率在10^‑3s量级，任意射向下的位置计算误差均小于5米，且计算结果在惯性系中表示，可直接参与弹上制导计算而无需额外的坐标转换。

Description

考虑高阶扰动引力影响的弹道偏差解析预报算法

技术领域

本发明涉及飞行动力学技术领域，特别地，涉及一种考虑高阶扰动引力影响的弹道偏差解析预报算法。

背景技术

自由飞行段是弹道导弹整个飞行段中飞行时间最长的一段，约占总飞行时间的90％以上。由于飞行高度较高，弹道导弹在自由飞行段主要受地球中心引力的作用，其飞行轨迹近似为椭圆轨道的一部分，但由于摄动因素(如地球非球型引力、稀薄大气阻力等)的存在，其真实轨迹又会偏离标准椭圆轨道。为保证弹道导弹的命中精度，在进行发动机关机控制时需要对摄动条件下的自由飞行段弹道偏差进行快速预报。实际上，考虑摄动因素的大气层外飞行器轨道预报是轨道动力学领域的经典问题之一，即初值问题。针对这一问题的经典理论主要有：平根数法、fg级数分解法、中间轨道法、非正交分解法以及自适应变步长数值积分方法等，但平根数法、fg级数分解法、中间轨道法、非正交分解法主要用于处理考虑J₂、J₃、J₄等低阶球谐项摄动的轨道预报问题，对高阶地球非球型引力摄动下的轨道预报问题则无法求解。自适应变步长数值积分法虽然可以计算高阶扰动引力影响下的轨道预报问题，但计算效率随着球谐阶次的增加迅速降低，无法应用于弹上实时计算。

因此，设计一种新的弹道偏差解析预报算法具有重要意义。

发明内容

本发明的目的在于提供一种考虑高阶扰动引力影响的弹道偏差解析预报算法，基于状态空间摄动法和沿飞行弹道的扰动引力重构模型提出。该方法的解算效率在10^-3s量级，任意射向下的位置计算误差均小于5米，且计算结果在惯性系中表示，可直接参与弹上制导计算而无需额外的坐标转换，具体技术方案如下：

一种考虑高阶扰动引力影响的弹道偏差解析预报算法，高阶扰动引力矢量分解和状态偏差解析预报模型推导；

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

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

其中，α_i均为常矢量系数，i＝0,1,…,9；表示真近点角为f时对应的标准二体弹道地心距，即且p表示二体弹道的半通径，e表示二体轨道偏心率；f为点P′在制导二体弹道上对应的真近点角，P′为制导二体弹道上与P点对应的点；

再获得扰动引力三分量的表达式为表达式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₅；

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

最后得到表达式6)：

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

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

式中：△v_r(f)、△v_β(f)和△v_z(f)分别为弹道状态偏差速度矢量在轨道柱坐标系中沿r轴、β轴和z轴方向的分量；△r(f)和△z(f)分别为弹道状态偏差位置矢量在轨道柱坐标系中沿r轴和z轴方向的分量；△t(f)为实际飞行时间与标准二体弹道飞行时间之差；h为二体弹道平面对应的动量矩矢量的模；为真近点角为ξ时对应的标准二体弹道地心距，即p为二体弹道的半通径，e为二体弹道的偏心率；

将表达式6)代入表达式7)中积分得到考虑扰动引力影响的自由段弹道偏差解析解如表达式8)-13)：

式中：ε₁-ε₄为

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)函数的解析表示式如下：

l^2,0,1(E)＝-a²(1-e²)cosE；

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

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

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

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

其中：插值核函数满足式中，ξ_k、和η_k为第k个节点的位置坐标，而ξ′_k、和η′_k分别为与该节点相邻的三个节点对应坐标轴方向上的坐标值，δg为弹道上任意一点的扰动引力矢量，δg_k为有限元网格八个节点的扰动引力矢量：

弹道上任意点P与轨道坐标系之间的几何关系图中，诸元二体弹道是基于标准关机点K_f确定的二体弹道，用于在导弹发射前进行扰动引力重构模型的构建；制导二体弹道是由导弹实际关机点K确定的二体弹道；实际弹道为导弹在地球非球型引力等摄动力作用下的真实飞行弹道；ψ₀和φ₀为地心角；P′为制导二体弹道上与P点对应的点，设P′在系中的位置矢量为则有表达式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在制导二体弹道上对应的真近点角；

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

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

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

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

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

附图说明

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

图1是实施例中弹道任意点位置坐标转换示意图；

图2(a)是实施例中由数值方法计算出的X方向位置偏差；

图2(b)是实施例中由解析方法计算出的X方向位置偏差与数值结果之差；

图3(a)是实施例中由数值方法计算出的Y方向位置偏差；

图3(b)是实施例中由解析方法计算出的Y方向位置偏差与数值结果之差；

图4(a)是实施例中由数值方法计算出的Z方向位置偏差；

图4(b)是实施例中由解析方法计算出的Z方向位置偏差与数值结果之差；

图5是不同条件下RKF数值方法与本实施例方法计算效率对比图。

具体实施方式

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

实施例：

一种考虑高阶扰动引力影响的弹道偏差解析预报算法，包括高阶扰动引力矢量分解和状态偏差解析预报模型推导，详情如下：

1、任意阶扰动引力矢量的分解，详情如下：

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

其中：插值核函数满足式中，ξ_k、和η_k为第k个节点的位置坐标，而ξ′_k、和η′_k分别为与该节点相邻的三个节点对应坐标轴方向上的坐标值，δg为弹道上任意一点的扰动引力矢量，δg_k为有限元网格八个节点的扰动引力矢量。

图1所示为弹道上任意点P与轨道坐标系之间的几何关系图，其中诸元二体弹道是基于标准关机点K_f确定的二体弹道，用于在导弹发射前进行扰动引力重构模型的构建；制导二体弹道是由导弹实际关机点K确定的二体弹道；实际弹道为导弹在地球非球型引力等摄动力作用下的真实飞行弹道；ψ₀和φ₀为地心角；P′为制导二体弹道上与P点对应的点，由于导弹在自由飞行段受到的摄动力为极小量，因此P′与P会非常接近。求解导弹自由飞行段弹道上任意点P的扰动引力矢量，需要首先计算该点在坐标系中的位置分量。在一阶假设下，P点处的扰动引力矢量需要基于P′点的位置进行计算。设P′在系中的位置矢量为则根据图1所示的几何关系则有表达式2)：

式中：为坐标系原点对应的地心距；Μ_z(·)和Μ_y(·)分别表示绕z轴和y轴变换的方向余弦矩阵，且有 表示真近点角为f时对应的标准二体弹道地心距，即且p表示二体弹道的半通径，e表示二体弹道的偏心率；△β、ψ₀和φ₀均为矢量夹角，△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在制导二体弹道上对应的真近点角；

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

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

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

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

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

定义变量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₅；

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

则有表达式6)：

2、状态偏差解析预报模型推导，具体是：

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

式中△v_r(f)、△v_β(f)和△v_z(f)分别为弹道状态偏差速度矢量在轨道柱坐标系中沿r轴、β轴和z方向的分量；△r(f)和△z(f)分别为弹道状态偏差位置矢量在轨道柱坐标系中沿r轴和z方向的分量；△t(f)为实际飞行时间与标准二体弹道飞行时间之差；h为二体弹道平面对应的动量矩矢量的模；表示真近点角为ξ时对应的标准二体弹道地心距，即p表示二体弹道的半通径，e表示二体弹道的偏心率；

将λ_j,s(f,ξ),j＝1,2,…,6,s＝1,3,5以及表达式6)代入表达式7)，并积分即获得考虑扰动引力影响的自由段弹道偏差解析解如表达式8)-13)：

式中：ε₁-ε₄为

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)函数的解析表示式如下：

l^2,0,1(E)＝-a²(1-e²)cosE；

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

其中：a表示二体轨道半长轴，e表示二体轨道偏心率，E表示偏近点角。

仿真实例：

通过数值仿真来验证本发明提出的弹道偏差解析预报模型的精度。假设地心惯性系中标准关机点位置矢量为x₀＝[0，6578140，0]^T，沿地心矢方向的初始速度为3300m/s，在导弹飞行平面内且垂直于地心矢方向上的初始速度为6680m/s。同时对不同射向进行遍历。导弹实际飞行中的关机点与标准关机点必定存在偏差，因此为验证本章提出方法的实际应用价值，仿真中对关机点状态进行拉偏处理，为此设置如表1所示的三组关机点状态偏差参数。

分别采用数值积分法、本节推导的解析解对弹道偏差进行解析预报，并将数值积分法的计算结果作为评价解析解精度的基准。仿真中扰动引力采用72×72阶EGM2008模型进行计算。

表1关机点状态偏差参数

序号	△x<sub>k</sub>(km)	△y<sub>k</sub>(km)	△z<sub>k</sub>(km)	△v<sub>xk</sub>(m/s)	△v<sub>yk</sub>(m/s)	△v<sub>zk</sub>(m/s)
							无拉偏	0	0	0	0	0	0
正向拉偏	100	100	100	50	50	10
							负向拉偏	-100	-100	-100	-50	-50	-10

仿真结果如图2-4所示，ρ_x，ρ_y和ρ_z分别为轨道坐标系中三个坐标轴方向上扰动引力引起的弹道偏差，该偏差值是由数值积分方法计算得到的，且扰动引力采用72×72阶球谐函数进行计算。△ρ_x，△ρ_y和△ρ_z分别为由本发明提出的解析方法计算出的弹道偏差与数值结果之差。结果表明在扰动引力引起的自由段弹道偏差超过700米的情况下，本发明提出的解析算法的位置偏差求解残差最大不超过2米，计算精度较高。

对本发明推导的扰动引力摄动下的弹道偏差解析预报模型的计算效率进行分析。多组仿真结果表明，提出的解析解进行一次弹道偏差预报需要耗时约3.8毫秒。为了体现该方法的效率，同样将其与RKF数值方法进行对比(见图5)。表2示为考虑72×72阶扰动引力条件下，RKF方法的积分精度和计算耗时随给定的绝对误差/相对误差直接的关系。可以看出，当给定的绝对/相对误差越小，积分精度约高，而积分耗时也越长。

表2考虑72×72阶扰动引力条件下RKF方法的积分精度和效率

绝对误差/相对误差	计算误差(m)	计算耗时(s)
			1.0×10<sup>-12</sup>	3.23960924×10<sup>-06</sup>	8.6700
1.0×10<sup>-11</sup>	2.94568389×10<sup>-05</sup>	5.7000
			1.0×10<sup>-10</sup>	2.90009027×10<sup>-04</sup>	3.6400
1.0×10<sup>-09</sup>	2.82104790×10<sup>-03</sup>	2.3450
			1.0×10<sup>-08</sup>	2.59484834×10<sup>-02</sup>	1.5650
1.0×10<sup>-07</sup>	0.21757124	1.1400
			1.0×10<sup>-06</sup>	1.44041035	0.7050
1.0×10<sup>-05</sup>	6.18215280	0.5900
			1.0×10<sup>-04</sup>	50.0562257	0.3900

图5所示为不同仿真条件下RKF数值积分方法与本文方法在计算效率方面的对比情况。横坐标为方法的位置积分误差，纵坐标为积分计算耗时。“五角星”表示本文提出的解析解的最大位置误差和计算耗时。可以看出：(1)在同等计算精度的情况下，本文提出的解析解的计算耗时为0.0038秒，且其计算时间与所考虑的扰动引力模型阶次无关；(2)不同扰动引力模型下RKF数值法的计算耗时不同，在同等计算精度条件下，8×8阶模型计算耗时约0.0235秒，36×36阶模型计算耗时约0.375秒，72×72阶模型计算耗时约1.14秒，108×108阶模型计算耗时约3.663秒，分别为本文提出解析解的6.2倍、98倍、300倍和963.9倍。显然，考虑的扰动引力阶次越高，解析解的相对计算效率越显著。

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

Claims (2)

1.一种考虑高阶扰动引力影响的弹道偏差解析预报算法，其特征在于：高阶扰动引力矢量分解和状态偏差解析预报模型推导；

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

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

其中，α_i均为常矢量系数，i＝0,1,…,9；表示真近点角为f时对应的标准二体弹道地心距，即且p表示二体弹道的半通径，e表示二体轨道偏心率；f为点P′在制导二体弹道上对应的真近点角，P′为制导二体弹道上与P点对应的点；

再获得扰动引力三分量的表达式为表达式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₅；

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

最后得到表达式6)：

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

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

式中：△v_r(f)、△v_β(f)和△v_z(f)分别为弹道状态偏差速度矢量在轨道柱坐标系中沿r轴、β轴和z轴方向的分量；△r(f)和△z(f)分别为弹道状态偏差位置矢量在轨道柱坐标系中沿r轴和z轴方向的分量；△t(f)为实际飞行时间与标准二体弹道飞行时间之差；h为二体弹道平面对应的动量矩矢量的模；表示真近点角为ξ时对应的标准二体弹道地心距，即p表示二体弹道的半通径，e表示二体弹道的偏心率；

λ_1,1(f,ξ)＝cos(f-ξ)；

λ_3,1(f,ξ)＝-sin(f-ξ)；

将表达式6)代入表达式7)中积分得到考虑扰动引力影响的自由段弹道偏差解析解如表达式8)-13)：

式中：ε₁-ε₄为

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

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

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

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

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

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

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

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

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

2.根据权利要求1所述的考虑高阶扰动引力影响的弹道偏差解析预报算法，其特征在于：P点处扰动引力矢量δg关于真近点角的函数表达式4)获得的具体过程如下：

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

其中：插值核函数满足式中，ξ_k、和η_k为第k个节点的位置坐标，而ξ′_k、和η′_k分别为与该节点相邻的三个节点对应坐标轴方向上的坐标值，δg为弹道上任意一点的扰动引力矢量，δg_k为有限元网格八个节点的扰动引力矢量：

弹道上任意点P与轨道坐标系之间的几何关系图中，诸元二体弹道是基于标准关机点K_f确定的二体弹道，用于在导弹发射前进行扰动引力重构模型的构建；制导二体弹道是由导弹实际关机点K确定的二体弹道；实际弹道为导弹在地球非球型引力等摄动力作用下的真实飞行弹道；ψ₀和φ₀为地心角；P′为制导二体弹道上与P点对应的点，设P′在系中的位置矢量为则有表达式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在制导二体弹道上对应的真近点角；

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

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

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

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