CN107766967B - 一种基于多项式拟合法的拦截弹预测制导方法 - Google Patents

Abstract

本发明涉及一种基于多项式拟合法的拦截弹预测制导方法，具体步骤如下：1、对拦截弹建模，列出拦截弹运动方程；2、计算拦截弹标准弹道，并列出拦截弹位移‑时间表；3、利用位移‑时间表对预测命中点进行粗算；4、改变初始弹道倾角，得到标准弹道族，不断迭代修正初始弹道倾角，修正预测命中点；5、根据预测命中点的位置，提出预测命中点制导方法，并在拦截弹飞行过程中继续修正预测命中点，降低脱靶量，若脱靶量不符合要求，继续使用下一步方法；6、采用多项式拟合法预报弹道的方法，并利用此方法，对预测命中点进一步修正，再次降低脱靶量，提高制导精度。本发明提高了制导精度，适应拦截弹在助推段过程中，周边环境不断剧烈变化的情况。

Description

一种基于多项式拟合法的拦截弹预测制导方法

技术领域

本发明涉及一种基于多项式拟合法的拦截弹预测制导方法，属于导弹制导领域，尤其涉及针对拦截弹预测命中点的解算、修正及拦截弹飞往预测命中点的制导方法。

背景技术

从1967年，时任美国总统的约翰逊下令部署“哨兵”***，美国开始了导弹防御***的研发。到1983年3月，里根政府提出发展导弹防御武器***的“战略防御倡议”(SDI)，该计划后来被称作“星球大战计划”。

今天，美国已经形成了包括陆基中段防御***(GBI)、海基中段防御***(SM-3)、末段高层防御***(THAAD)和近程防御***(PAC-3)在内的全空域导弹防御***。这些防御***均采用动能杀伤战斗部(KV)，要求制导***具有非常高的精度。其中，GBI、SM-3和THAAD均为中远程地空导弹。在进入末制导之前，导弹需要飞行相当长的时间。传统的比例导引法在这个阶段将不再适用，需要采用基于预测命中点的制导方法。

C.Grubin在1964年提出了用于火箭跟踪实现设计好的弹道的算法，里面涉及了预测制导的基本思想。M.Salama在1987提出了预测制导的概念，主要针对目标运动的预测，并且比较了预测制导与纯追踪法的拦截弹道和需要过载。对目标弹道的预测是预测制导方法的一个重要组成部分。X.Zhang总结了前期弹道导弹弹道预报的方法。随着临近空间飞行器概念的发展，新的目标特性引起研究的兴趣。L.Qin研究了临近空间飞行器HTV-2的弹道预报方法。另一方面，拦截弹自身弹道的预报也是研究的一个重要方向。如G.CHatterji研究了短期弹道的预报方法。L Hainz研究了拦截弹在线弹道预报的方法，随后提出了改进的线性化方法用于快速弹道预报。随着制导律的研究进展,一些高级的制导律需要用到从当前时刻到拦截时刻之间所谓“剩余飞行时间”。关于“剩余飞行时间”出现了很多研究成果。C.Tournes将预测制导方法与二阶滑模控制算法结合起来。拦截弹弹道预报方面的最新进展是BT Burchett使用高斯伪谱法进行快速弹道预报。

发明内容

本发明的目的是针对传统的比例导引法或者增强比例导引法容易产生过大的法向过载，影响最终的拦截效果，提出一种基于多项式拟合法的拦截弹预测制导方法。

本发明为一种基于多项式拟合法的拦截弹预测制导方法，具体步骤如下：

1、对拦截弹建模，列出拦截弹运动方程；

2、计算拦截弹标准弹道，并列出拦截弹位移-时间表；

3、利用位移-时间表对预测命中点进行粗算；

4、改变初始弹道倾角，得到标准弹道族，不断迭代修正初始弹道倾角，修正预测命中点；

5、根据预测命中点的位置，提出预测命中点制导方法，并在拦截弹飞行过程中继续修正预测命中点，降低脱靶量，若脱靶量不符合要求，继续使用下一步方法；

6、采用多项式拟合法预报弹道的方法，并利用此方法，对预测命中点进一步修正，再次降低脱靶量，提高制导精度。

本发明与现有技术相比，具有的有益效果是：提高了制导精度，适应拦截弹在助推段过程中，周边环境不断剧烈变化的情况。

附图说明

图1是本发明运行流程图。

图2是零攻角飞行的标准弹道。

图3是大气层内交战场景。

图4是标准弹道族，从左至右初始弹道倾角为80°至30°(每隔10°绘制一条弹道)。

图5是拦截弹初始弹道倾角选择。

图6是单次预测命中点计算的弹道倾角。

图7是单次预测命中点计算的拦截弹道。

图8是单次预测命中点计算的弹道倾角。

图9是单次预测命中点计算的拦截弹道。

图10是偏差弹道。

图11是标准弹道上的弹道倾角。

图12是标准弹道上的速度。

图13是速度拟合结果。

图14是弹道倾角拟合结果。

图15是纵程拟合结果。

图16是高度拟合结果。

图17是标准弹道与偏差弹道。

具体实施方式

下面结合附图和实例，对本发明的技术方案做进一步的说明。

本发明是一种基于多项式拟合法的拦截弹预测制导方法，如图1所示，具体包括如下步骤：

1、对拦截弹建模，列出拦截弹运动方程

二维平面中，拦截弹零攻角飞行的动力学方程如下：

其中，v为拦截弹飞行速度，g为重力加速度，γ为俯仰角，t为飞行时间，ρ₁为空气密度，C_D为阻力系数，S为拦截弹参考面积，m为拦截弹质量，x为水平方向飞行距离，y为竖直方向飞行高度。

2、计算拦截弹标准弹道，并列出拦截弹位移-时间表

在没有外界干扰因素的情况下，若拦截弹严格按照零攻角飞行，则拦截弹的飞行状态与初始状态以及飞行时间t是一一对应的关系。每一个初始状态对应一条不同的弹道。本发明中，将每一个初始状态对应的零攻角弹道定义为标准弹道。另外，将每个时刻t拦截弹的位置与拦截弹的初始位置之间的距离定义为位移，用R(t)表示。

图2给出了一条零攻角飞行的标准弹道。该标准弹道初始条件为：高度(15KM)、速度(2000m/s)、弹道倾角(45deg)。表1给出了该标准弹道上不同飞行时刻对应的位移。3、利用位移-时间表对预测命中点进行粗算

考虑如图3所示交战场景：

首先给定一个飞行时间t_n，目标从B点飞行到达P点。为了击中目标，拦截弹需要经过时间t_n从A点飞行到达P点。

令P点到A点的距离为r，拦截弹经过时间t_n能够飞行的位移为R(t_n)。其中，R(t_n)可以通过表1初始弹道倾角45°的位移-时间表插值得到。若R(t_n)＜r，说明拦截弹无法经过时间t_n到达P点，需要增加飞行时间t_n；反之，则需要减少飞行时间t_n。若R(t_n)＝r，则P点即为预测命中点。

表1

下面计算一个算例：令拦截弹初始位置为(0，15km)，初始速度为2000m/s。目标初始位置为(120km，60km)，初始速度为2000m/s，水平向左做匀速直线运动。通过迭代，可以得到预测命中点P的位置为(37.42km，60km)。令预测命中点P到拦截弹初始位置A点的连线与水平面的夹角为预测命中点的高低角θ_P，如图3中所示，

θ_P与初始弹道倾角不同。我们可以根据新的θ_P计算新的位移-时间表，通过新的位移-时间表，可以插值得到R_50.2(41.2891)＝59.079km，与前面用45°初始弹道倾角标准弹道的位移-时间表插值得到的R₄₅(41.2891)＝58.523km之间有较大的误差，从而使得预测命中点计算产生较大的误差。为了解决这个问题，在预测命中点计算中，不能使用单一的位移-时间表。

4、利用标准弹道族，修正预测命中点

不同的初始弹道倾角对应的标准弹道上，相同的飞行时间对应位移有较大的差别。为了满足不同高低角的预测命中点的计算需求，需要在数据库中存储多条不同初始弹道倾角的标准弹道。本发明将这样的一组标准弹道定义为“标准弹道族”。

图4给出了一组不同初始弹道倾角的标准弹道。图中标准弹道的初始弹道倾角间隔为10°。在实际计算中，存储在数据库的标准弹道的初始弹道倾角间隔为0.1°。

使用初始弹道倾角45°的位移-时间表计算得到的预测命中点为P₁(37.42km，60km)，该预测命中点的高低角

为了提高预测命中点的计算精度，使用初始弹道倾角50.2°的位移-时间表重新计算预测命中点。

可以得到预测命中点P₂的位置为(37.852km，60km)，飞行时间t_go＝41.07s。

通过计算，可以得到预测命中点P₂的高低角为

同样地，再次使用初始弹道倾角49.9°的位移-时间表重新计算预测命中点，可以得到预测命中点P₃的位置为(37.832km，60km)，飞行时间t_go＝41.08s。

通过计算，可以得到预测命中点P₃的高低角为

因为数据库中存储的标准弹道的初始弹道倾角间隔为0.1°，所以到此迭代结束。

5、介绍预测命中点制导方法，并在拦截弹飞行过程中继续修正预测命中点

首先，前面预测命中点计算结果中，预测命中点的高低角θ_P并不能作为拦截弹的初始弹道倾角。若拦截弹以预测命中点的高低角θ_P作为初始弹道倾角进行无攻角飞行，拦截弹并不能到达预测命中点，如图5所示。因此需要根据预测命中点的高低角θ_P来确定拦截弹的初始弹道倾角。这仍然是一个迭代的过程。这样，就确定了拦截弹的初始弹道倾角。拦截弹以此初始弹道倾角做无攻角飞行，即可到达预测命中点。

然而在实际飞行过程中，拦截弹从当前弹道倾角变到指令弹道倾角需要一个过程。因此，实际飞行过程中，应该将标准弹道上的实时的弹道倾角作为指令。

图6给出了单次预测命中点计算的指令弹道倾角与实际飞行弹道的弹道倾角。其中指令弹道倾角为标准弹道上相应时刻的弹道倾角。图7给出了单次预测命中点计算的拦截弹道，脱靶量为1308.3m。

为了提高拦截精度，可以采用多次计算命中点的方法，比如每隔1s计算一次预测命中点，直到t_go＜10s。

图8给出了多次预测命中点计算的指令弹道倾角与实际飞行弹道的弹道倾角。其中指令弹道倾角为标准弹道上相应时刻的弹道倾角。图9给出了多次预测命中点计算的拦截弹道，脱靶量为119.5m。

6、介绍利用多项式拟合法预报弹道的方法，并利用此方法，对预测命中点进一步修正

只有极少部分的微分方程具有解析解，大部分微分方程初值问题可以通过数值方法如龙格库塔法等进行求解。在某些需要对快速性要求较高的场合，可以使用拉格朗日多项式拟合的方法近似求解微分方程。

标准弹道初始时刻的状态参数为x₀，y₀，v₀，γ₀，t时刻的状态参数为x₁(t)，y₁(t)，v₁(t)，γ₁(t)。现在，拦截弹在初始时刻的状态参数相对标准弹道的状态参数有一个变换量Δx，Δy，Δv，Δγ，如图10所示。下面求解拦截弹t时刻的状态参数x₂(t)，y₂(t)，v₂(t)，γ₂(t)。

标准弹道上的状态参数满足如下微分方程：

同样，真实弹道上的状态参数也满足微分方程

设真实弹道与标准弹道之间的状态参数偏差为δx(t)、δy(t)、δv(t)和δγ(t)

δv(t)＝v₂(t)-v₁(t) (18)

δγ(t)＝γ₂(t)-γ₁(t) (19)

δx(t)＝x₂(t)-x₁(t) (20)

δy(t)＝y₂(t)-y₁(t) (21)

将真实弹道上的状态参数的微分方程组与标准弹道上的状态参数的微分方程组分别相减，并忽略二阶以上小量，可以得到

上述微分方程组为关于状态变量δx、δy、δv和δγ的时变线性微分方程组，时变参数x₁(t)、y₁(t)、v₁(t)、γ₁(t)为标准弹道上的状态变量，它们的值存储在标准弹道数据库里。该微分方程组可以用多项式拟合近似求解。由于微分方程组中，

和

与状态变量δx无关，因此不需要第3个微分方程就能将状态变量δv、δγ和δy求解出来。

将状态变量δv、δγ和δy用三阶多项式拟合，如下

δv(t)＝a₁t³+b₁t²+c₁t+δv₀ (26)

δγ(t)＝a₂t³+b₂t²+c₂t+δγ₀ (27)

δy(t)＝a₃t³+b₃t²+c₃t+δy₀ (28)

对上式两边进行求导，可以得到

选择时刻t₁、t₂和t₃，满足微分方程(22)、(23)、(25)，即

上述方程为九元一次方程组，确定时刻t₁、t₂和t₃以及标准弹道上对应时刻的状态参数v₁(t₁)、v₁(t₂)、v₁(t₃)、γ₁(t₁)、γ₁(t₂)、γ₁(t₃)、y₁(t₁)、y₁(t₂)和y₁(t₃)，就可以求出拟合多项式的参数a₁、b₁、c₁、a₂、b₂、c₂、a₃、b₃和c₃。

下面，以一个实例来说明具体的计算流程。

这里，我们将一个初始速度v₀＝2000，初始弹道倾角γ₀＝45^°，初始位置为(0，15000)的惯性飞行弹道作为标准弹道。

这条标准弹道上的速度v₁和弹道倾角λ₁随时间的变化如图11和图12所示。在总的飞行时间60s以内，任何时刻的状态参数弹道倾角γ₁和速度v₁都可以通过插值得到。

假设拦截弹初始位置和初始弹道倾角不变，初始速度变为v＝990，即δv＝-10。接着，我们分步求解新的速度和弹道倾角。

1)确定时间节点：t₁＝20，t₂＝40，t₃＝60；

2)插值求解时间节点上的弹道倾角γ₁、速度v₁和高度y₁：v₁(20)＝1372.7、v₁(40)＝1233.4、v₁(60)＝1136.8、γ₁(20)＝0.6914、γ₁(40)＝0.5694、γ₁(60)＝0.4238、y₁(20)＝3.6072×10⁴、y₁(40)＝5.1421×10⁴和y₁(t₃)＝6.2741×10⁴；

3)解(32)至(40)方程组，求得拟合多项式的系数：a₁＝1.0993×10^-6，b₁＝3.0669×10^-4，c₁＝-0.0225，a₂＝9.6198×10^-10，b₂＝-5.6703×10^-7，c₂＝-2.2894×10^-5，a₃＝3.0518×10^-5，b₃＝-0.0055，c₃＝-7.1137。

于是，速度、弹道倾角和高度变化量为

δv(t)＝1.0993×10^-6t³+3.0669×10^-4t²-0.0225t-10 (41)

δγ(t)＝9.6198×10^-10t³-5.6703×10^-7t²-2.2894×10^-5t (42)

δy(t)＝3.0518×10^-5t³-0.0055t²-7.1137t (43)

在此基础上，就可以求解新的弹道。与前面步骤类似，将状态变量δx和δy用三阶多项式拟合，如下

δx(t)＝a₄t³+b₄t²+c₄t+δx₀ (44)

对上式两边进行求导，可以得到

在时刻t₁、t₂和t₃，满足微分方程(24)即

接着，我们分步求解新的弹道。

1)确定时间节点：t₁＝20，t₂＝40，t₃＝60；

2)插值求解时间节点上的弹道倾角γ₁、速度v₁和高度y₁：v₁(20)＝1372.7、v₁(40)＝1233.4、v₁(60)＝1136.8、γ₁(20)＝0.6914、γ₁(40)＝0.5694、γ₁(60)＝0.4238、y₁(20)＝3.6072×10⁴、y₁(40)＝5.1421×10⁴和y₁(t₃)＝6.2741×10⁴；

根据δv(t)、δγ(t)和δy(t)的拟合多项式，求解弹道倾角和速度以及高度的偏差：δv(20)＝-9.9674、δv(40)＝-9.4846、δv(60)＝-8.3641，δγ(20)＝-0.0014、δγ(40)＝-0.0043、δγ(60)＝-0.0078，δy(20)＝-144.2、δy(40)＝-291.3、δy(60)＝-439.9；

3)解方程组(46)至(48)

求得拟合多项式的系数：a₄＝3.8776×10^-5，b₄＝-0.0080，c₄＝-7.0832。

于是，纵程变化量为

δx(t)＝3.8776×10^-5t³-0.0080t²-7.0832t (49)

图13至图16为速度、弹道倾角、纵程、和高度的拟合结果，可见拟合的结果较好。

实际飞行过程中，拦截弹的状态参数与标准弹道上相同时刻的状态参数不一致，因此预测命中点还需要进一步修正。下面以一个具体的算例来说明预测命中点修正的方法。

拦截弹初始位置为(0，15000)，初始速度v_M0＝1990m/s。目标初始位置为(135000，60000)，初始速度v_T0＝2000m/s。

根据之前的方法计算出的预测命中点为(44960，60000)，预测命中点的高低角为44.97°，标准弹道的初始弹道倾角为50.4°，剩余飞行时间t_go＝45.02s。若拦截弹的初始速度v_M0＝2000m/s，与标准弹道上的初始速度相同，最后的脱靶量为MD＝8.31m；若拦截弹的初始速度与标准弹道上的初始速度有偏差，为v_M0＝1990m/s，则脱靶量为MD＝193.28m。两次拦截的弹道如图17所示。

为了进一步提高预测命中点的精度，采用多项式拟合的方法对预测命中点进行修正。

根据上一部分的方法，可以计算出拦截弹的初始速度为v_M0＝1990m/s时，位置的偏差为

δx(t)＝4.5048×10^-5t³-0.0077t²-6.4308t (50)

δy(t)＝4.0991×10^-5t³-0.0067t²-7.8348t (51)

根据位置偏差方程，就可以得到新的位移-时间表，并重新计算预测命中点。

下面给出计算流程。

1)利用50.4°标准弹道以及位置偏差方程计算得到新的预测命中点P₁(44603,60000)，预测命中点高低角

通过插值计算得到经过时间t_go之后，拦截弹的位置与初始位置连线的高低角为θ₁＝44.94°；

2)将标准弹道更换为初始弹道倾角为θ_v2＝θ_v1+θ_P1-θ₁＝50.7°的标准弹道。根据新的标准弹道计算新的位置偏差方程。

δx₂(t)＝4.5064×10^-5t³-0.0077t²-6.3911t (52)

δy₂(t)＝4.1277×10^-5t³-0.0067t²-7.8718t (53)

根据新的标准弹道和新的位置偏差方程，重新计算预测命中点P₂(44623,60000)，预测命中点高低角

通过插值计算得到经过时间t_go之后，拦截弹的位置与初始位置连线的高低角为θ₂＝45.24°。

与θ₂偏差足够小，停止迭代。

通过上面的计算流程，得到初始速度v_M0＝1990m/s时，预测命中点P(44623,60000)，标准弹道50.7°。最终脱靶量为MD＝36.22m。

Claims (1)

1.一种基于多项式拟合法的拦截弹预测制导方法，其特征在于：该方法具体步骤如下：

步骤一、对拦截弹建模，列出拦截弹运动方程；

步骤二、计算拦截弹标准弹道，并列出拦截弹位移-时间表；

步骤三、利用位移-时间表对预测命中点进行粗算；

步骤四、改变初始弹道倾角，得到标准弹道族，不断迭代修正初始弹道倾角，修正预测命中点；

步骤五、根据预测命中点的位置，提出预测命中点制导方法，并在拦截弹飞行过程中继续修正预测命中点，降低脱靶量，若脱靶量不符合要求，继续使用下一步方法；

步骤六、采用多项式拟合法预报弹道，并利用此方法，对预测命中点进一步修正，再次降低脱靶量，提高制导精度；具体过程如下：使用拉格朗日多项式拟合的方法近似求解微分方程；

标准弹道初始时刻的状态参数为x₀，y₀，v₀，γ₀，t时刻的状态参数为x₁(t)，y₁(t)，v₁(t)，γ₁(t)；现在，拦截弹在初始时刻的状态参数相对标准弹道的状态参数有一个变换量Δx，Δy，Δv，Δγ，下面求解拦截弹t时刻的状态参数x₂(t)，y₂(t)，v₂(t)，γ₂(t)；

标准弹道上的状态参数满足如下微分方程：

同样，真实弹道上的状态参数也满足微分方程

设真实弹道与标准弹道之间的状态参数偏差为δx(t)、δy(t)、δv(t)和δγ(t)

δv(t)＝v₂(t)-v₁(t) (18)

δγ(t)＝γ₂(t)-γ₁(t) (19)

δx(t)＝x₂(t)-x₁(t) (20)

δy(t)＝y₂(t)-y₁(t) (21)

将真实弹道上的状态参数的微分方程组与标准弹道上的状态参数的微分方程组分别相减，并忽略二阶以上小量，可以得到

上述微分方程组为关于状态变量δx、δy、δv和δγ的时变线性微分方程组，时变参数x₁(t)、y₁(t)、v₁(t)、γ₁(t)为标准弹道上的状态变量，它们的值存储在标准弹道数据库里；该微分方程组可以用多项式拟合近似求解；由于微分方程组中，

和

与状态变量δx无关，因此不需要第3个微分方程就能将状态变量δv、δγ和δy求解出来；

将状态变量δv、δγ和δy用三阶多项式拟合，如下

δv(t)＝a₁t³+b₁t²+c₁t+δv₀ (26)

δγ(t)＝a₂t³+b₂t²+c₂t+δγ₀ (27)

δy(t)＝a₃t³+b₃t²+c₃t+δy₀ (28)

对上式两边进行求导，可以得到

选择时刻t₁、t₂和t₃，满足微分方程(22)、(23)、(25)，即

上述方程为九元一次方程组，确定时刻t₁、t₂和t₃以及标准弹道上对应时刻的状态参数v₁(t₁)、v₁(t₂)、v₁(t₃)、γ₁(t₁)、γ₁(t₂)、γ₁(t₃)、y₁(t₁)、y₁(t₂)和y₁(t₃)，就可以求出拟合多项式的参数a₁、b₁、c₁、a₂、b₂、c₂、a₃、b₃和c₃。

Images