一种考虑导航误差的最优拦截制导方法

1.本发明属于飞行器制导控制技术领域，具体涉及考虑导航误差的最优拦截制导方法。


背景技术：

2.自1957年世界上第一颗人造卫星发射成功后，人类便真正地进入了太空时代。而后的几十年间，随着航天技术不断发展，越来越多的卫星被送入太空，并且在通信、导航、气象观测以及科学研究等领域发挥了重要的作用，给人类的生活带来了极大的便利。与此同时，各类卫星在军事领域的应用也在快速发展。在现代战争中，由于卫星的轨道高度较高，不受国界和地理条件限制，且不易受到攻击，因此其成为获取战场信息的主要来源。为了在对抗中建立优势，有必要研究如何对这些卫星进行有效拦截。
3.针对拦截起点固定而拦截终点不固定的单脉冲能量最优拦截问题，传统的方法往往是针对目标器和拦截器的标称轨道进行设计，进而在任务约束时间范围内采用一维搜索方法获得能量最优的拦截时刻。但由于实际任务中导航误差的存在，若按照标称设计轨道进行拦截，则终端拦截精度很有可能无法满足末制导的初始条件要求，从而导致拦截任务以失败告终。因此，在拦截过程中，必须通过施加相应的修正脉冲去减小导航误差带来的影响。


技术实现要素：

4.本发明的目的是为了解决现有实际任务中导航误差的存在，若按照标称设计轨道进行拦截，则终端拦截精度很有可能无法满足末制导的初始条件要求，从而导致拦截任务以失败告终的问题，而提出一种考虑导航误差的最优拦截制导方法。
5.一种考虑导航误差的最优拦截制导方法具体过程为：
6.步骤一、给定初始时刻t0以及初始时刻拦截器和目标器的标称轨道参数后，设定t0为拦截任务的起始时刻，在任务给定的拦截时间范围内采用一维搜索算法寻找能量最优的终端拦截时刻并计算得到初始时刻施加的脉冲
7.步骤二、以步骤一中得到的能量最优解作为参考，在初始时刻导航误差存在的情况下，基于最优拦截条件求解最优拦截时刻和初始时刻施加脉冲的误差，同时给出相应的终端拦截误差；
8.步骤三、在给定修正脉冲的施加时刻t1后，采用解析方法估计施加修正脉冲的幅值以及修正后的终端拦截误差，从而确定修正脉冲幅值以及修正后的终端拦截误差与修正时刻t1的解析关系；
9.步骤四、确定修正脉冲施加时刻或对应的时间范围。
10.本发明的有益效果：
11.本发明设计了一种考虑导航误差的最优拦截制导方法，以弥补以往只针对标称轨
道参数设计工作的不足。该方法考虑了目标器和拦截器的导航误差对能量最优拦截过程的影响，并解析给出了修正脉冲的估计方法，进而设计了一种脉冲修正策略，该策略有利于拦截任务的顺利实现。
12.本发明提出了一种考虑导航误差的最优拦截制导方法。在本发明中，基于最优拦截条件和相关的偏导数矩阵，解析地给出了最优拦截时刻和初始时刻施加脉冲的误差与初始导航误差的关系。进一步地，对于给定的修正脉冲施加时刻，给出了修正脉冲幅值和终端拦截误差的解析估计方法。采用本发明方法求解导航误差下的能量最优拦截问题，只需要选择合适的修正脉冲时刻，便可以满足任务要求的终端拦截精度。此外，本发明也适用于需要多次施加修正脉冲的情况，只需要在执行完成当前的修正脉冲后再采用本发明中的方法计算下一次满足任务要求的修正脉冲即可。
附图说明
13.图1为考虑导航误差的lambert最优制导算法流程图；
14.图2为单脉冲能量最优拦截问题的几何示意图；
15.图3为修正脉冲幅值和终端拦截误差标准差随修正时刻的变化关系图；
16.图4为加权指标随修正时刻的变化关系图。
具体实施方式
17.具体实施方式一：本实施方式一种考虑导航误差的最优拦截制导方法具体过程为：
18.在本方法中所采用的动力学模型均为二体模型，表示为：
[0019][0020][0021]
其中，μ表示地球的引力常数，r和v分别表示j2000惯性系航天器的位置矢量和速度矢量，|r|则表示相应位置矢量的大小；表示r的一阶导数，表示v的一阶导数；
[0022]
步骤一、给定初始时刻t0以及初始时刻拦截器和目标器的标称轨道参数后，设定t0为拦截任务的起始时刻，在任务给定的拦截时间范围内采用一维搜索算法寻找能量最优的终端拦截时刻并计算得到初始时刻施加的脉冲
[0023]
步骤二、以步骤一中得到的能量最优解作为参考(拦截时刻以及相应的脉冲这两者是对应的，已知已知便可通过lambert算法求得)，考虑在初始时刻导航误差存在的情况下，基于最优拦截条件求解最优拦截时刻和初始时刻施加脉冲的误差(公式13和15)，同时给出相应的终端拦截误差(公式16)；
[0024]
步骤三、在给定修正脉冲的施加时刻t1后，采用解析方法估计施加修正脉冲的幅值(以脉冲的幅值标准差评价)以及修正后的终端拦截误差(用终端拦截误差标准差衡量)，从而确定修正脉冲幅值以及修正后的终端拦截误差与修正时刻t1的解析关系(公式24和25即为相应的关系)；
[0025]
步骤四、按照设定的某种指标最优或者根据任务给定的相关精度要求确定合适的修正脉冲施加时刻或对应的时间范围。
[0026]
具体实施方式二：本实施方式与具体实施方式一不同的是，所述步骤一中给定初始时刻t0以及初始时刻拦截器和目标器的标称轨道参数后，设定t0为拦截任务的起始时刻，在任务给定的拦截时间范围内采用一维搜索算法寻找能量最优的终端拦截时刻并计算得到初始时刻施加的脉冲
[0027]
具体过程为：
[0028]
所述拦截器可以是在轨运行的己方拦截卫星，也可以是己方发射入轨后的拦截导弹；
[0029]
所述目标器是敌方在轨运行的卫星或其他航天器；
[0030]
给定初始时刻t0和拦截时间范围[t
fmin
,t
fmax
]，初始时刻拦截器和目标器的标称轨道参数在j2000地心惯性坐标系下给出，初始时刻拦截器和目标器的标称轨道参数包括拦截器初始时刻的位置矢量和速度矢量以及目标器初始时刻的位置矢量和速度矢量；
[0031]
其中拦截器初始时刻的位置矢量和速度矢量分别记为r
10
和v
10
，目标器初始时刻的位置矢量和速度矢量分别记为r
20
和v
20
；
[0032]
当给定拦截时间范围内的任意一个拦截时刻tf后，则目标器终端位置矢量r2可通过求解kepler方程得到。记初始位置矢量r1＝r
10
，则转移轨道所需的初始速度矢量v
t1
便可通过求解lambert问题得到，即v
t1
＝f(r1,r2,δt)，其中δt＝t
f-t0。也就是说，初始时刻所施加的脉冲幅值大小是拦截时间的一元函数，即δv1＝g(tf)。因此，初始时刻施加的最小的脉冲便可通过一维搜索拦截时间tf得到。需要说明的是，由于初始时刻是预先给定的，可假设t0＝0，则得δt＝tf，因此在本发明中认为二者含义一致，不进行区分。
[0033]
在任务给定的拦截时间范围内采用一维搜索算法寻找能量最优的终端拦截时刻并计算得到初始时刻施加的脉冲
[0034]
其它步骤及参数与具体实施方式一相同。
[0035]
具体实施方式三：本实施方式与具体实施方式一或二不同的是，所述在任务给定的拦截时间范围内采用一维搜索算法寻找能量最优的终端拦截时刻并计算得到初始时刻施加的脉冲具体过程为：
[0036]
在步骤一中，采用的一维搜索方法是分段黄金分割法，具体求解过程描述如下：
[0037]
首先将拦截时间范围[t
fmin
,t
fmax
]等分为n个子区间，在每个子区间内采用黄金分割法求解对应的能量最优解
[0038]
最后，通过寻找所得的n个能量最优解中的最小值即可获得整个拦截时间范围内的能量最优解(这两者是对应的，求出便可通过lambert算法求得)；
[0039]
式中，为第k个子区间内的最优拦截时刻，为第k个子区间内对应的初始时刻施加的脉冲，为拦截时间范围[t
fmin
,t
fmax
]内能量最优的终端拦截时刻，为拦截时间范围[t
fmin
,t
fmax
]内能量最优的初始时刻施加的脉冲。
[0040]
其它步骤及参数与具体实施方式一或二相同。
[0041]
具体实施方式四：本实施方式与具体实施方式一至三之一不同的是，所述步骤二中以步骤一中得到的能量最优解作为参考，考虑在初始时刻导航误差存在的情况下，基于最优拦截条件求解最优拦截时刻和初始时刻施加脉冲的误差(公式13和15)，同时给出相应的终端拦截误差；
[0042]
具体过程为：
[0043]
在步骤二中，最优拦截条件可表示为
[0044][0045]
其中
[0046][0047]
式中，δv1为初始时刻所需要施加的速度脉冲，δt为最优拦截时间，v
t1
为转移轨道所需的初始速度矢量，t为求转置；
[0048]
同时，目标器终端位置矢量(拦截位置)r2可表示为
[0049]
r2＝f
tr20
g
tv20
ꢀꢀ
(3)
[0050]
其中
[0051][0052]
式中,p2表示目标器轨道的半通径，和分别表示目标器初始时刻和拦截时刻的真近点角，μ表示地球的引力常数；f
t
和g
t
为拉格朗日系数；|r2|表示目标器终端位置矢量的大小；
[0053]
对式(3)关于拦截时间δt求偏导，则可得
[0054][0055]
最终可得到目标器终端位置矢量(拦截位置)r2关于拦截时间δt的偏导数为
[0056][0057]
其中
[0058][0059]
式中，i3表示3阶单位矩阵；
[0060]
目标器拦截时刻真近点角关于拦截时间δt的偏导数为
[0061][0062]
式中，a2和e2分别为目标轨道的半长轴和偏心率，e2为目标器在拦截时刻的偏近点角；
[0063]
假设拦截器在初始时刻的位置和速度导航误差为[δr
10
,δv
10
]，目标器在初始时刻的位置和速度导航误差为[δr
20
,δv
20
]；
[0064]
则在导航误差存在的情况下，初始时刻所需要施加的速度脉冲为
[0065][0066]
式中，r1表示拦截器初始位置矢量，δ(δt)表示最优拦截时间的误差，δr2表示目标器终端位置矢量误差；
[0067]
是通过求解lambert问题的偏导数矩阵得到的(参见zhang g,zhou d,mortari d,et al.covariance analysis of lambert's problem via lagrange's transfer-time formulation[j].aerospace science and technology,2018,77:765-773.)(已知r1、r2、δt，将r1、r2、δt作为lambert问题的输入，得到v
t1
；进而在已知r1、r2、δt、v
t1
的情况下，采用上述论文中的方法可输出)；
[0068]
δr2则可表示为
[0069][0070]
式中，φ
(t)
表示目标器从初始时刻到拦截时刻的二体状态转移矩阵，其求解方法已比较成熟(一种简单的方法可参见reynolds r g.direct solution of the keplerian state transition matrix[j].journal of guidance control and dynamics,2022.doi:10.2514/1.g006373.)；
[0071]
同时，该二体状态转移矩阵还可表示为分块矩阵的形式
[0072]
[0073]
式中，表示目标器拦截时刻位置矢量对初始时刻位置矢量的偏导数，表示目标器拦截时刻的位置矢量对初始时刻速度矢量的偏导数，表示目标器拦截时刻的速度矢量对初始时刻位置矢量的偏导数，表示目标器拦截时刻的速度矢量对初始时刻速度矢量的偏导数；φ
(t)
为6
×
6的矩阵，均表示3
×
3的矩阵；
[0074]
根据最优拦截条件，可得
[0075][0076]
于是，可求得最优拦截时间的误差为
[0077][0078]
其中
[0079][0080]
式中，k表示转移轨道所需的初始速度矢量v
t1
对最优拦截时间δt的全导数；
[0081]
为了方便描述，记
[0082][0083]
式中，a1表示δ(δt)受δr
10
影响的系数矩阵，a2表示δ(δt)受δv
10
影响的系数矩阵，a3表示δ(δt)受δr
20
影响的系数矩阵，a4表示δ(δt)受δv
20
影响的系数矩阵；
[0084]
则最优拦截时间的误差可重新写为
[0085]
δ(δt)＝a1δr
10
a2δv
10
a3δr
20
a4δv
20
ꢀꢀ
(13)
[0086]
初始时刻施加脉冲(速度增量)的误差为
[0087][0088]
同样，可简记为
[0089]
δ(δv1)＝b1δr
10
b2δv
10
b3δr
20
b4δv
20
ꢀꢀ
(15)
[0090]
其中
[0091][0092]
式中，b1表示δ(δv1)受δr
10
影响的系数矩阵，b2表示δ(δv1)受δv
10
影响的系数矩阵，b3表示δ(δv1)受δr
20
影响的系数矩阵，b4表示δ(δv1)受δv
20
影响的系数矩阵；
[0093]
最终可得到在导航误差存在的情况下，拦截时刻拦截器和目标器的相对位置误差为
[0094][0095]
其中
[0096][0097]
式中表示拦截器在拦截时刻的位置矢量，表示拦截器在拦截时刻的位置误差矢量，c1表示δ(δrf)受δr
10
影响的系数矩阵，c2表示δ(δrf)受δv
10
影响的系数矩阵，c3表示δ(δrf)受δr
20
影响的系数矩阵，c4表示δ(δrf)受δv
20
影响的系数矩阵；
[0098]
需要说明的是的求解过程与的求解过程基本一致，只需要将相应的参数更换为转移轨道的参数即可(公式3-8是基于目标器的参数，即r
20
，v
20
和δt进行推导的，其他参数都可在已知这三个量之后计算得到，而转移轨道对应的这三个参数为r
10
(或r1，前面定义了二者是等价的)，v
t1
和δt，将r
20
，v
20
和δt相应替换为r
10
，v
t1
和δt，再进行公式3-8的推导便可求解得到)；
[0099]
由此，若给定初始时刻拦截器的导航误差协方差矩阵为[p
r10
,p
v10
]，目标器的导航误差协方差矩阵为[p
r20
,p
v20
]，则最优拦截时刻的误差方差为：
[0100][0101]
最优拦截时刻的误差标准差为：
[0102][0103]
初始时刻施加脉冲的误差协方差矩阵为
[0104][0105]
初始时刻施加脉冲幅值的标准差为
[0106][0107]
式中，trace()表示求解矩阵的迹；
[0108]
终端相对位置的误差协方差矩阵为
[0109][0110]
终端相对距离的标准差为
[0111]
[0112]
式中，e[
·
]表示数学期望；
[0113]
在步骤二和步骤三中，终端拦截误差具体为拦截器和目标器在最优拦截时刻的相对距离误差。
[0114]
其它步骤及参数与具体实施方式一至三之一相同。
[0115]
具体实施方式五：本实施方式与具体实施方式一至四之一不同的是，所述步骤三中在给定修正脉冲的施加时刻t1后，采用解析方法估计施加修正脉冲的幅值(以脉冲的幅值标准差评价)以及修正后的终端拦截误差(用终端拦截误差标准差衡量)，从而确定修正脉冲幅值以及修正后的终端拦截误差与修正时刻t1的解析关系(公式24和25即为相应的关系)；具体过程为：
[0116]
假设修正时刻得到的拦截器和目标器的位置速度导航误差分别为[δr
11
,δv
11
]和[δr
21
,δv
21
]；
[0117]
此时，按照修正时刻得到的拦截器和目标器的位置速度导航误差分别为[δr
11
,δv
11
]和[δr
21
,δv
21
]可预测得到目标器在终端拦截时刻的位置误差为
[0118][0119]
式中，为目标器从修正时刻t1到拦截时刻tf的状态转移矩阵；为目标器拦截时刻位置矢量对修正时刻位置矢量的偏导数矩阵，为目标器拦截时刻位置矢量对修正时刻速度矢量的偏导数矩阵；
[0120]
假设t1时刻施加的修正脉冲矢量为q，则按照修正时刻得到的拦截器和目标器的位置速度导航误差分别为[δr
11
,δv
11
]和[δr
21
,δv
21
]计算得到拦截器在终端拦截时刻的位置误差为
[0121][0122]
式中，为拦截器从修正时刻t1到拦截时刻tf的状态转移矩阵；为拦截器拦截时刻位置矢量对修正时刻位置矢量的偏导数矩阵，为拦截器拦截时刻位置矢量对修正时刻速度矢量的偏导数矩阵；
[0123]
由于施加修正脉冲的目的是使得拦截器和目标器在拦截时刻处于同一位置，即满足如下条件
[0124][0125]
则可计算得到施加的修正脉冲为
[0126][0127]
上式可以简记为
[0128]
[0129]
其中
[0130][0131]
式中，是d1表示q受δr
10
影响的系数矩阵；d2表示q受δv
10
影响的系数矩阵；d3表示q受δr
20
影响的系数矩阵；d4表示q受δv
20
影响的系数矩阵；d5表示q受δr
11
影响的系数矩阵；d6表示q受δv
11
影响的系数矩阵；d7表示q受δr
21
影响的系数矩阵；d8表示q受δv
21
影响的系数矩阵；
[0132]
则施加修正脉冲的协方差矩阵为
[0133]
pq＝e[qq
τ
]
[0134]
施加修正脉冲幅值的标准差为
[0135][0136]
在此基础上，拦截器的真实初始状态，在施加通过带有导航误差的轨道数据计算得到的“不准确”的初始脉冲和修正脉冲后，在最优拦截时刻，拦截器和目标器的真实的相对位置误差可表示为
[0137][0138]
式中，是拦截器真实的初始状态在施加“不准确”的初始脉冲和修正脉冲后，在拦截时刻的位置误差；
[0139]
则施加修正后的终端相对位置误差协方差矩阵为
[0140][0141]
施加修正后的终端相对距离的标准差为
[0142][0143]
综上估计得到在t1时刻施加的修正脉冲q和施加修正后的终端相对位置误差以及计算得到相应修正脉冲幅值标准差σ(q)和修正后的终端拦截误差标准差
[0144]
本发明中，采用终端拦截误差的1σ值作为衡量终端拦截精度的指标，其中σ(δrf)对应于不施加修正脉冲的情况，对应于施加修正脉冲的情况；
[0145]
其它步骤及参数与具体实施方式一至四之一相同。
[0146]
具体实施方式六：本实施方式与具体实施方式一至五之一不同的是，所述步骤四中按照设定的某种指标最优或者根据任务给定的相关精度要求，确定合适的修正脉冲施加
时刻或对应的时间范围；具体过程为：
[0147]
在步骤四中，某种指标可假设为终端拦截误差标准差和修正脉冲幅值标准差σ(q)的加权指标。根据步骤三中的结果可知和σ(q)均是修正脉冲施加时刻t1的一元函数，可表示为
[0148][0149]
式中，h1()表示终端拦截误差标准差和修正脉冲施加时刻的函数关系，h2()表示修正脉冲幅值标准差和修正脉冲施加时刻的函数关系；
[0150]
因此，可设计加权指标j为
[0151][0152]
其中j表示加权指标，w1和w2分别表示两个指标的权重系数，且满足w1 w2＝1，和是归一化参数，可根据任务需求合理设计使得归一化后两个指标的量级相当；
[0153]
最终，加权指标可表示为修正脉冲施加时刻的一元函数，即j＝h3(t1)，进而可通过一维搜索算法确定使得加权指标j最优的修正脉冲施加时刻；
[0154]
式中，h3()表示加权指标j和修正脉冲施加时刻t1的函数关系。
[0155]
其它步骤及参数与具体实施方式一至五之一相同。
[0156]
具体实施方式七：本实施方式与具体实施方式一至五之一不同的是，所述步骤四中按照设定的某种指标最优或者根据任务给定的相关精度要求，确定合适的修正脉冲施加时刻或对应的时间范围；具体过程为：
[0157]
在步骤四中，给定任务的相关精度要求可以考虑为对的精度要求，如设置的精度上限为
[0158]
由于是脉冲修正时刻的一元函数。因此，可通过二分法求解方程的零根，从而得到满足终端拦截精度要求的修正脉冲施加时刻范围；
[0159]
式中，h1()表示终端拦截误差标准差和修正脉冲施加时刻的函数关系。
[0160]
其它步骤及参数与具体实施方式一至五之一相同。
[0161]
采用以下实施例验证本发明的有益效果：
[0162]
实施例一：
[0163]
设初始时刻拦截器初始时刻的标称位置速度矢量为：
[0164][0165]
目标器初始时刻的标称位置速度矢量为：
[0166][0167]
拦截器的导航误差标准差为
[0168][0169]
目标器的导航误差标准差为
[0170][0171]
首先，设置初始时刻t0＝0，拦截终端时刻的搜索范围为[0,14400]s，然后在标称状态下，求解拦截起点固定而拦截终点不定的能量最优单脉冲拦截问题。最终通过分段黄金分割法搜索得到最优拦截时刻标称值为初始时刻施加脉冲的标称值为其脉冲幅值大小为
[0172]
在考虑初始时刻拦截器和目标器的导航误差存在的情况下，根据步骤二中的解析方法可计算得到最优拦截时刻、机动脉冲幅值以及终端拦截误差标准差，并采用monte-carlo方法打靶仿真1000次进行验证，结果如表1所示：
[0173]
表1初始导航误差对最优拦截参数的影响(1σ)
[0174][0175]
从上述结果可以发现，在考虑初始导航误差存在的情况下，若按标称设计结果进行拦截，终端拦截误差标准差(1σ)约为50km，拦截任务极有可能失败。因此，必须在拦截过程中施加相应的修正脉冲弥补导航误差对终端拦截精度的影响。
[0176]
首先根据步骤三中的结果，分析不同修正脉冲施加时刻对施加脉冲幅值和终端拦截误差的影响，具体情况如图3所示；
[0177]
情况一：设置归一化参数加权系数取为w1＝w2＝0.5。则加权指标随修正时刻的变化情况如图4所示；
[0178]
图4显示，加权指标在整个拦截过程中存在一个极小值点，因此采用黄金分割法求得加权指标极小值处对应的修正时刻为t1＝5333.767s。采用步骤三中的方法可估计得到该时刻施加的修正脉冲幅值标准差为σ(q)＝0.0142km/s，修正后的终端拦截误差标准差为为了验证所得结果的正确性，采用monte-carlo方法打靶仿真1000次，对所得结果进行统计分析，得到修正脉冲幅值标准差为0.0144km/s，修正后的终端拦截误差标准差为13.7425km。可以看出与解析方法所得结果基本一致，证明了本发明中方法的有效性。
[0179]
情况二：设置任务对终端拦截精度的要求为由于
为拦截时刻t1的函数，因此，只要求解方程的零根即可。采用二分法求解得到方程的解为t1＝6332.827s，即修正脉冲只需在该时刻之后添加，即可满足任务给定的终端拦截精度要求。由本发明中的解析方法可计算得到，在该时刻施加的修正脉冲幅值标准差为σ(q)＝0.0200km/s，修正后的终端拦截误差标准差为采用monte-carlo方法打靶仿真1000次，对所得结果进行统计分析，得到修正脉冲幅值标准差为0.0203km/s，修正后的终端拦截误差标准差为10.0493km。由此说明，本发明所提出的方法是正确有效的。
[0180]
本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。