为应对高超声速导弹的威胁，美国洛·马公司很早就开始了对THAAD系统的改进方案的研究——THAAD增程型(THAAD-ER)。THAAD增程型助推器由一级变为两级，极大地提高了助推时间和关机速度。2018年7月19日，洛·马公司又签下价值1.64亿美元的合同，负责升级在韩国部署的THAAD系统。升级后的THAAD系统射程可能增加至3倍，拦截区域增加至9~12倍，即使初始拦截失败，也可以有足够时间继续进行拦截，提高成功率。
自德国V-2导弹研制成功以来，对导弹的拦截一直是各国重视的难题，各类制导律也应运而生。比例导引法是使用最为广泛的导引方法之一，诸多研究人员也以比例导引法为基础进行了大量研究。Zarchan^[1]对比例导引法的推导过程和线化形式进行了详细分析。Kumar等^[2]对比例导引法进行了改进，考虑拦截弹在助推段速度急剧上升，在滑翔段速度也有较大变化，用线化模型求解零控脱靶量和剩余飞行时间会产生很大误差，因此在求解时选用非线性模型，再通过数值积分进行计算，同时通过最优控制的方法推导出比例导引法的增益随时间变化的函数，将比例导引法中恒定的增益改为可变的增益，并根据拦截弹飞行阶段的不同给出不同的表达式，有效提高了比例导引法的性能。李辕等^[3-4]使用线化的比例导引法，对顺轨拦截和逆轨拦截2种模式分别给出了剩余飞行时间的估计方法。Dhananjay^[5]和Ghosh^[6]等在假设导弹使用比例导引法飞行的情况下，给出了剩余飞行时间的估计方法。
由于THAAD增程型拦截弹本身结构和特点的限制，在进入末制导之前导弹有较长时间无控滑翔段，中制导完全依赖于助推段，传统制导律如比例导引法需要在拦截过程的全程进行控制以应对目标的机动，无法满足拦截的需求，预测制导方法则因能够很好解决这一问题而得到广泛应用。Grubin等^[7]在1964年提出了跟踪设计好的弹道算法，制导方式利用存储好位置、速度的文件及剩余飞行时间对火箭进行控制，使火箭在主动段结束时到达入轨点和入轨速度，涉及了预测制导的基本思想。Salama等^[8]提出了预测制导的概念，要求导弹的速度向量时刻指向预测命中点，比较了预测制导与纯追踪法的拦截弹道和需要过载，突出了预测制导的优势。张华伟等^[9]通过选取拦截高度，解出目标由飞行至拦截高度所需时间及位置坐标，并计算出拦截弹飞行至拦截点所需时间，不断迭代解算出预测命中点。Zhang等^[10]用迭代飞行时间的方法计算预测命中点，但估计拦截弹的剩余飞行时间时仅简单地利用弹目斜距除以接近速度。Zhao等^[11]用较为简单的模型进行计算预测命中点，假设拦截弹运动至预测命中点的弹道接近直线。Hahn等^[12]在高空对拦截弹的运动方程进行线化，得到拦截弹弹道的解析式，进而求解拦截时间及脱靶量，并施加控制消除脱靶量。Dwivedi等^[13]研究目标的广义状态估计和模型预测制导，制导律以对零控脱靶量的预测为基础，同时对零控脱靶量和剩余飞行时间进行预测，使结果更加精确，但需要进行数值积分运算。Song^[14]和Zhang^[15]等利用神经网络方法训练得到制导指令，拦截弹弹道接近最优弹道，取得较好的拦截效果，但若要做到在拦截弹射程内进行全范围拦截，则需要进行大量的神经网络训练工作。陈万春^[16]和尤刘球^[17]在标准弹道中插值求得剩余飞行时间，解算预测命中点完成拦截任务。
文献[10-12]在计算预测命中点时使用的弹目相对运动模型非常简单，在线化后很容易得出剩余飞行时间和预测命中点。但当考虑的模型较为复杂(如飞行距离较长)时，拦截弹弹道不再近似为直线或者拦截弹速度变化较大则会产生很大的误差。文献[13]在计算剩余飞行时间和预测命中点时较为准确，但需要进行大量的数值积分运算，会对弹载计算机造成沉重的负担，并且当建立的拦截弹模型更贴近实际、更为复杂时，会增大误差。文献[14-15]使用训练神经网络的方法获得剩余飞行时间及指导指令，但一次的训练结果仅能满足一个小范围内的拦截任务，若想完成在拦截弹射程内的全范围拦截，需要进行大量的训练工作，难于实现。文献[16-17]在计算预测命中点时也基于拦截弹弹道近似直线的假设，当拦截射程较远时，该方法无法使用，并且需在弹上装载大量的标准弹道文件，对弹上存储空间的要求较高。
因此，需要一种既能适应复杂的拦截弹模型又能完成远距离拦截任务的预测制导方法。本文基于标准弹道族的思想，提出了一种通过解析方法计算剩余飞行时间进而求解预测命中点的方法，避免了迭代过程，计算准确，运算效率高。
1 拦截弹建模 根据网上公开数据进行测绘反设计，绘制如图 1所示的拦截弹外形图。拦截弹总长6.2 m，一级助推器直径0.53 m，二级助推器直径0.37 m。

图 1 拦截弹主要尺寸 Fig. 1 Main dimensions of interceptor

图选项

拦截弹发射质量为1 100 kg，发射采用70°俯仰角倾斜发射。弹体由两级固体助推器和动能拦截器(KKV)组成，一二级助推器均可进行推力矢量控制，助推器参数如表 1所示，末制导阶段释放KKV完成动能拦截，KKV的导引头安装在侧面，可通过直接力控制KKV的姿态完成制导^[18]。本文研究拦截弹中制导阶段。
表 1 助推器参数 Table 1 Parameters of boosters

级数	质量/kg	推力/kN	燃烧时间/s
1	860	80	30
2	170	30	8

表选项






在地面坐标系下对拦截弹建立动力学和运动学模型，其原点一般位于制导雷达所在位置，x轴位于水平面内并指向正北方向，y轴垂直地面指向天空方向，z轴按右手法则决定。
2 基于标准弹道的预测制导方法 本文主要研究拦截弹的预测制导问题，不对目标的弹道预报进行分析，目标预报弹道视为已知。如图 2所示，拦截任务的实质是使拦截弹与目标同时到达同一点。本文的预测制导方法核心在于：先求解剩余飞行时间t_go，求得t_go后目标到达的点即预测命中点，之后求得预测命中点所在标准弹道并控制拦截弹按照标准弹道飞行，完成拦截任务。标准弹道是在给定固定制导指令后，在线下通过数值积分得到的拦截弹弹道，因此无论拦截弹模型多么复杂，只要在拦截仿真时给拦截弹以同样的指令就可以使拦截弹弹道与标准弹道相同。使用标准弹道的预测制导流程如图 3所示^[19-21]。

图 2 拦截弹对抗目标示意图 Fig. 2 Schematic diagram of interceptor against target

图选项

图 3 预测制导流程[19-21] Fig. 3 Predictive guidance flowchart[19-21]

图选项

已知目标弹道后，通过在目标弹道上利用时间插值可得到预测命中点，因此只有准确计算出拦截弹飞行时间，才能得到准确的预测命中点。本文以标准弹道为基础，通过迭代飞行时间法^[16]和迭代预测命中点法2种方法分别求解预测命中点，并对2种方法进行对比。
标准弹道是本文的核心工具，下面首先介绍标准弹道。
2.1 标准弹道的生成 考虑二维平面内一二级连续助推的情况。拦截弹的动力学和运动学方程可写为^[22]

(1)

(2)

(3)

(4)

(5)

(6)

(7)

式中：m为质量；V为速度；P为推力；α为攻角；X为在空气中受到的阻力；G为重力；θ为弹道倾角；J为转动惯量；ω为弹体转动角速度；M为所受力矩；为俯仰角；x和y分别为地面坐标系下的纵程和高度。
主动段结束后，拦截弹处于无控飞行状态，其状态与飞行时间一一对应，被动段的初始状态是决定弹道的唯一因素，本文将这样一条弹道定义为标准弹道，拦截弹位置与发射点之间的连线定义为拦截弹飞行的斜距，用R表示。
2.2 标准弹道族 为使标准弹道能覆盖拦截范围内的空间点，需要生成一簇标准弹道，称这样一簇标准弹道为标准弹道族。在生成标准弹道时，助推器连续助推，不考虑一二级级间点火间隔。
如图 4所示，所有从(0, 0)出发的均为标准弹道，本文中将标准弹道进行编号，标准弹道的编号为被动段初始弹道倾角的100倍，用N_θ表示。图 4中标准弹道间隔100号，在实际计算中，数据库中相邻标准弹道的间隔为1号即0.01°。与标准弹道相交的为等时间线，等时间线从40 s开始，每相邻2条时间线时间间隔为20 s，直至220 s结束。得到标准弹道族和等时间线后，当确定预定的拦截位置时，可以在图中通过测绘手段找到拦截弹所需跟踪的标准弹道及飞行时间。但在实际拦截任务中，需要在很短时间内完成这一工作，人工测绘显然不可行，需要寻求标准弹道族和飞行时间的数学规律，进而形成算法编程实现。

图 4 标准弹道族 Fig. 4 A bunch of standard trajectories

图选项

2.3 迭代飞行时间法 图 5为一个简单拦截场景。图中：AM为标准弹道，BT为目标弹道目标从B向T飞行，A点为拦截弹发射点。首先考虑拦截弹在初始发射时刻的情况，此时剩余飞行时间t_go与飞行时间t_f相等，计算步骤如下^[16]：

图 5 拦截示意图 Fig. 5 Schematic diagram of interception

图选项

步骤1?? 设定初始飞行时间t_f，任选一条标准弹道N_θ，拦截弹按照标准弹道飞行，目标按照已知的预报弹道飞行。
步骤2?? 如果经过t_f后目标到达P点、拦截弹到达M点，设拦截弹距离发射点的斜距为R_M(t_f)，目标距离拦截弹发射点的斜距为R_T(t_f)，则有以下3种情况：
1) 如果经过t_f后目标到达P′、拦截弹到达M′，R_M(t_f)与R_T(t_f)相等，t_f为所求的飞行时间，P′点即为预测命中点，程序结束。
2) 如果经过t_f后目标到达P₁点、拦截弹到达M₁点，R_M(t_f) < R_T(t_f)，表明t_f过短，需增加t_f，重复进行步骤2。
3) 如果经过t_f后目标到达P₂点、拦截弹到达M₂点，R_M(t_f)>R_T(t_f)，表明t_f过长，需减少t_f，重复进行步骤2。
上述是发射时刻的特例，此时t_f与剩余飞行时间t_go相同，在拦截弹飞行过程中往往还需继续计算t_go，在计算R_M(t_f)与R_T(t_f)时，只需将起点由发射点A改为拦截弹当前位置，所得即为剩余飞行时间，之后在目标弹道上插值求得预测命中点。
由于选取的标准弹道为任选的一条标准弹道，首次计算的预测命中点往往不在标准弹道上，导弹按照初始选取的标准弹道飞行无法到达预测命中点以完成拦截任务。因此，需要选取其他适当的标准弹道，重复上面的运算，使得到的预测命中点既在目标弹道也在标准弹道上，这样拦截弹按照标准弹道飞行即能完成拦截任务。
2.4 选取标准弹道 在文献[16]中已经对这一问题提出了较好的解决方法，但规划的标准弹道近似直线，预测命中点的解算是基于弹道近似直线的假设，在拦截距离较近的拦截弹(如THAAD拦截弹)上效果很好，但无法应用于远距离拦截任务中，在进行远程拦截时需要寻找新的方法。下面用一个简单的例子说明。
THAAD拦截弹的最大拦截距离200 km，弹道近似直线，如图 6所示。预测命中点与发射点间连线AP₁与水平面的夹角θ₁和被动段初始弹道倾角θ₀近似相等，因此，可通过求解θ₁求得粗算的预测命中点所在的标准弹道来进行标准弹道的迭代。设预测命中点的坐标为(x_P, y_P)，拦截弹发射点的坐标为(x_A, y_A)，则所求标准弹道的代号N_θ可表示为

图 6 不同射程拦截对比 Fig. 6 Comparison of interception with different ranges

图选项

(8)

THAAD增程型拦截弹拦截距离为THAAD拦截弹的3倍左右，最大射程可达600 km，若采用文献[16]的标准弹道计算策略，则会有较大误差。
当拦截点坐标为(400, 120)km时，如图 7所示，标准弹道被动段初始弹道倾角为

图 7 远距离拦截示例 Fig. 7 An example of long-range interception

图选项

(9)

若拦截弹以此弹道倾角为指令，在射程100 km附近时，拦截弹弹道仍近似直线弹道。随着飞行距离和飞行时间增长，由于飞行高度低拦截弹速度损失严重，高度明显下降，拦截弹无法到达预定拦截位置。
为解决这一问题，在拦截距离较远时选用高抛弹道。利用高抛弹道后，就带来了新的问题，如图 6所示，AP₂与地平面的夹角θ₃与弹道倾角θ₂没有明显关联，无法通过现有的方法得到标准弹道代号选取新的标准弹道。
在主动段结束后，拦截弹由于只受空气动力和重力的影响，弹道近似于抛物线，若能通过求解预测命中点P₂所在的抛物线，进而就能确定标准弹道。
首先求解标准弹道对应的近似抛物线。抛物线方程可表示为

(10)

在标准弹道上选取3点，即可求得未知数a、b、c。求解出每条标准弹道对应的a、b、c值，可以得到一组a、b、c数据。将参数a、b、c表示为标准弹道的代号N_θ的函数，则抛物线方程就可以表示为

(11)

如果能求得参数a、b、c与N_θ间的函数关系，就可以通过将预测命中点P的坐标(x_P, y_P)代入式(11)求得对应的N_θ。
本文利用多项式拟合法通过最小二乘法原理求解^[23]。设多项式方程为以下形式：

(12)

由最小二乘法的定义可知，当已知k个点数据(x, y)，需要确定各系数d使得各点的偏差e的平方和最小：

(13)

把k个点数据代入式(12)，可得

(14)

将式(14)各系数用矩阵表示，令

(15)

(16)

(17)

则式(14)可表示为AD=Y，等式两端左乘A^T可得

(18)

求解此线性方程组，可得解为d′₀, d′₁, …, d′_n，拟合多项式为

(19)

若使用三阶多项式对a、b、c进行拟合，则拟合结果可表示为

(20)

式(11)就可以表示为

(21)

将预测命中点坐标(x_P, y_P)代入式(21)，可解出N_θ，即确定了标准弹道。这样就解决了已知预测命中点后求其所在的标准弹道的问题，再通过迭代飞行时间法即可解出满足精度要求的预测命中点。
由于有空气动力的影响，拦截弹的被动段弹道并不是严格意义上的抛物线，用一条抛物线拟合仅能满足部分区域拟合精度较高，可以根据需求进行分段拟合。本文按拦截弹射程分为0~200 km、200~400 km、400~600 km三段进行拟合，抛物线近似情况如图 8所示，图中的标准弹道为间隔200号绘制，拟合抛物线与标准弹道拟合情况很好。由于标准弹道是在地面坐标系中建立的，由于地表为球面，随着射程的增长，y坐标与实际高度的偏差逐渐增大，为能更直观地了解高度信息，在图 8中还绘制了地球表面。

图 8 抛物线近似标准弹道情况 Fig. 8 Parabolic approximation of standard trajectory

图选项

图 9~图 11为a、b、c三个系数不同射程下的拟合情况。可知，拟合情况很好。

图 9 系数a、b、c拟合情况(射程：0~200 km) Fig. 9 Fitting of coefficient a, b and c (Range: 0-200 km)

图选项

图 10 系数a、b、c拟合情况(射程：200~400 km) Fig. 10 Fitting of coefficient a, b and c (Range: 200-400 km)

图选项

图 11 系数a、b、c拟合情况(射程：400~600 km) Fig. 11 Fitting of coefficient a, b and c (Range: 400-600 km)

图选项

2.5 迭代预测命中点法 2.3节中已经通过迭代飞行时间法完成了解算预测命中点的问题，但该方法存在的明显问题给使用者带来了难度：①需要使用者通过经验给定标准弹道和剩余飞行时间的初值，算法对初值较为敏感，使用者需要通过长期的经验积累才能选取适当的初值；②需要在弹上装载大量弹道文件，并不断进行迭代插值运算，给弹载计算机造成了很大负担，程序运行时间较长。拦截过程中，在已知目标预报弹道的情况下，期望拦截位置是一个直观概念，如果根据期望的拦截位置设定初始条件，可以极大降低算法使用的难度。因此，设定期望的拦截位置为初始预测命中点，之后再迭代预测命中点来确定最终的预测命中点进行预测制导，这时就需要通过其他方法得到飞行时间。下面通过继续挖掘标准弹道数据，用解析方法计算剩余飞行时间，不仅可以使初值设定更容易，还可以减少大量的迭代和插值运算步骤，提高程序效率。
由2.4节可知，已知预测命中点(x_P, y_P)后，可以求得该点所在的标准弹道N_θ，问题就可以简化为：已知(x_P, y_P)和N_θ，求解对应的拦截弹飞行时间t_f，再减去当前时刻求得剩余飞行时间t_go。当拦截距离较远时，高抛弹道中拦截弹高度y随时间先上升后下降，不是单调函数，而射程x随时间单调增加，因此考虑使用多项式拟合法求解射程x与飞行时间t的函数关系。在x方向上，拦截弹在被动段只受空气阻力的作用，且在飞行高度较高时速度损失很小，速度大小、高度及空气密度变化小，x方向上空气阻力变化较小，因此选用二次多项式拟合x-t函数关系。

(22)

如果再拟合出a₂、b₂、c₂与N_θ的函数关系，就能通过N_θ和x求解飞行时间t_f。再次利用多项式拟合法，可拟合出a₂、b₂、c₂与N_θ的函数关系：

(23)

x-t函数关系可以表示为

(24)

解得

(25)

在根的取舍上，因为是实际物理问题，一定有且仅有唯一根满足条件。分析式(24)，由物理意义可知，a₂(N_θ) < 0。从数学关系考虑，随着t逐渐增加，x先上升后下降，因此有2个根满足数学关系。但实际情况中，x一定为单调上升的，因此仅有较小的一个根为有物理意义的根。

(26)

利用抛物线近似x-t情况如图 12所示。k、m、n分别取9、5、5，在N_θ=1 600、N_θ=3 000两处进行分段拟合，图 13为拟合情况，可知拟合情况很好。

图 12 抛物线近似x-t情况 Fig. 12 Parabolic approximation of x-t

图选项

图 13 系数a2、b2、c2拟合情况 Fig. 13 Fitting of coefficient a2, b2 and c2

图选项

性能配置为：CPU为Intel Core i7-4790、主频3.60 GHz，运行内存4.00 GB的64位计算机，以VC++6.0软件平台，在相同的模型下运行拦截程序。表 2为迭代预测命中点法与迭代飞行时间法运行时间对比。由于预测制导方法仅在拦截弹助推段进行，选取拦截弹助推段程序运行时间进行对比。
表 2 程序运行时间对比 Table 2 Program run time comparison

仿真次数	设定迭代飞行时间法运行时间/s	设定迭代飞行时间法脱靶量/m	设定预测命中点法运行时间/s	设定预测命中点法脱靶量/m	运行时间节省比例/%	脱靶量减少比例/%
1	28.72	96.89	23.21	16.21	19.19	83.27
2	29.17	293.83	23.29	83.85	20.16	71.46
3	29.22	79.10	23.39	79.10	19.95	0
4	29.14	296.38	22.90	174.68	21.41	41.06
5	29.32	99.86	22.45	79.49	23.43	20.40
6	29.03	178.72	23.26	199.82	19.88	-11.81
7	28.88	251.91	23.12	311.36	19.94	-23.60

表选项






表 2中每次仿真均取相同拦截位置利用2种不同方法进行拦截仿真。迭代预测命中点法在保证脱靶量的前提下，运行时间减少20%左右，能够有效提升程序运行效率，迭代预测命中点法优于迭代飞行时间法。下文的拦截仿真使用迭代预测命中点法进行。
3 拦截仿真结果分析 在VC++6.0平台建立六自由度模型，使用迭代预测命中点法进行拦截仿真，地面坐标系的原点选在目标弹道弹下点。拦截弹一二级助推器连续助推，选择迭代预测命中点法在助推段实时计算更新预测命中点，预测命中点更新间隔取1 s。当拦截弹与目标间距离达到最小时停止仿真，此距离为脱靶量^[24-25]。
3.1 不同射程 拦截弹发射阵地布置在目标弹下点上，在目标弹道上取距离拦截弹发射阵地不同水平距离的点作为初始设定的预测命中点。拦截弹脱靶量随射程增加有逐渐增大的趋势，在射程600 km时达到最大385.45 m，脱靶量多数情况下小于100 m。拦截弹终点位置与设定的预测命中点，平均偏差333 m，最大偏差1 161 m，相对于射程而言偏差很小，完成在预计拦截点附近进行拦截的要求。拦截仿真结果如表 3所示。
表 3 不同射程拦截仿真结果 Table 3 Simulation results of interception with different ranges

设定的预测命中点			射程/km	脱靶量/m	终点			拦截点偏移/m	末速度/ (m·s-1)
x坐标/m	y坐标/m	z坐标/m			x坐标/m	y坐标/m	z坐标/m
-34 862	48 202	-35 956	50.09	15.79	-34 870	48 188	-35 965	19	4 686.31
-69 770	47 692	-71 961	100.02	19.06	-69 626	47 677	-71 812	208	4 066.02
-110 274	46 608	-113 736	158.42	28.24	-110 274	46 636	-113 736	29	3 164.37
-138 417	45 542	-142 763	198.54	36.85	-138 200	45 514	-142 540	312	2 781.90
-173 091	43 878	-178 526	248.31	81.13	-172 846	43 810	-178 273	359	2 587.59
-218 956	41 080	-225 831	314.38	23.63	-218 838	41 112	-225 709	173	2 566.10
-243 346	39 316	-250 987	349.59	63.27	-243 350	39 379	-250 991	63	2 598.21
-277 268	36 542	-285 973	398.33	86.28	-277 275	36 628	-285 981	86	2 660.54
-320 014	32 516	-330 062	460.45	194.60	-320 514	32 270	-330 577	759	2 740.87
-349 067	29 441	-360 027	502.05	199.82	-349 473	29 596	-360 446	604	2 796.76
-386 483	25 077	-398 618	554.99	1.57	-386 325	25 095	-398 455	228	2 838.07
-432 395	19 101	-445 971	620.12	385.45	-431 661	19 588	-445 214	1 161	2 834.60

表选项






图 14为射程分别为600 km和50 km时的拦截仿真结果，网格代表地球表面。图 15为纵平面内拦截仿真情况，虚线为目标弹道，实线为不同射程下拦截弹弹道。图 16为拦截弹速度变化曲线，拦截弹末速度均大于2 500 m/s，保持较大动能。图 17为拦截弹在不同位置拦截同一目标时高度随时间变化的曲线，当射程较近时，拦截弹为近似直线弹道，在拦截弹上升段拦截目标，与THAAD拦截弹类似，当拦截距离较远时，为高抛弹道，利用拦截弹下降段拦截目标。

图 14 三维拦截仿真图例(射程：50 km及600 km) Fig. 14 Examples of 3D interception simulation (Range: 50 km and 600 km)

图选项

图 15 纵平面内拦截情况 Fig. 15 Interception in vertical plane

图选项

图 16 拦截弹速度 Fig. 16 Velocity of interceptor

图选项

图 17 拦截弹高度 Fig. 17 Altitude of interceptor

图选项

3.2 不同航路捷径 航路捷径指导弹发射点至空中目标的航向在水平面上投影的垂直距离，是拦截弹拦截能力的一个重要指标。航路捷径对脱靶量的影响对验证制导方法的可行性有着重要意义。
使用与上文相同的地面坐标系，即坐标原点依然为弹下点，设定4个不同的初始预测命中点，改变航路捷径进行仿真分析，点坐标如表 4所示，结果如表 5所示。脱靶量平均值为145.37 m，存在航路捷径时，预测制导方法依然有效。
表 4 设定的预测命中点坐标 Table 4 Expected predicted impact point coordinate

序号	x/m	y/m	z/m
1	-34 862	48 202	-35 956
2	-69 770	47 692	-71 961
3	-138 417	45 542	-142 763
4	-277 268	36 542	-285 973

表选项






表 5 存在航路捷径时的脱靶量 Table 5 Miss distance with course shortcut

航路捷径/km	脱靶量/m
	点1	点2	点3	点4
0	3.64	19.06	36.85	86.28
50	6.05	4.85	33.76	118.10
100	36.12	82.94	24.95	134.39
150	178.31	75.37	33.20	258.81
200	87.18	124.92	71.07	449.50
250	55.86	42.33	127.59	415.55
300	37.56	52.74	153.13	326.55
350	95.33	104.96	217.85	135.17
400	180.90	191.41	353.20	128.68
450	260.17	230.11	225.70
500	260.78	168.21	109.11
550	125.82	103.70	167.82
600	319.18	377.57

表选项






图 18为拦截仿真的弹下点轨迹，虚线为目标弹下点，实线为拦截弹弹下点。图中：由上至下分别对应表 4中4个设定的预测命中点，为保证图片清晰，绘图时仅取部分航路捷径情况的结果(选取的航路捷径见图中标尺)。

图 18 拦截弹与目标的弹下点 Fig. 18 Projective point of interceptor and target

图选项

3.3 气动不确定性对仿真结果的影响 拦截弹的气动系数一般可以通过气动拟合公式、DATCOM软件等简单得到，但用这些方法得到的气动系数必然会与真实情况存在一定偏差，分析预测制导方法因气动系数不准确而产生的脱靶量偏差十分必要。在气动系数中加入噪声进行100次仿真，噪声均值为气动系数大小的5%，仿真结果如表 6和图 19所示。
表 6 不确定性仿真结果 Table 6 Uncertainty simulation results

射程/km	无噪声脱靶量/m	有噪声平均脱靶量/m	有噪声脱靶量标准差/m	平均飞行时间/s	有噪声脱靶量最大值/m
158.42	28.24	27.94	2.08	69.73	32.66
349.59	63.27	83.21	11.25	147.65	118.10
554.99	1.57	18.04	13.02	210.51	60.86
620.12	385.45	384.89	23.31	228.39	435.42

表选项

图 19 不确定性仿真的脱靶量 Fig. 19 Miss distance of uncertainty simulation

图选项

加入噪声后的平均脱靶量相比无噪声时的脱靶量增加不多，100次打靶仿真中的脱靶量最大值相比无噪声时有一定增加。随着射程增加，飞行时间逐渐增加，由于拦截弹的特殊结构，助推段时间固定，若飞行时间增长，则无控飞行的时间增长，气动系数的变化产生的效应会产生累积，造成脱靶量的波动增大，与仿真结果中脱靶量标准差随着射程增加而增加相符。
图 20为使用迭代预测命中点法进行预测制导的程序运行流程。

图 20 预测制导程序流程 Fig. 20 Flowchart of predictive guidance program

图选项

4 结论 本文针对THAAD增程型拦截弹提出了基于标准弹道的预测制导方法，并通过仿真分析验证预测制导方法的能力，主要取得以下成果：
1) 根据网络公开数据，基于反导拦截的各阶段工作原理，对THAAD增程型拦截弹进行建模，并编写了六自由度拦截仿真程序。
2) 规划了高抛弹道扩充标准弹道族，完成大射程拦截任务，使用抛物线分段近似标准弹道，利用多项式拟合法计算抛物线参数，求得抛物线系数关于标准弹道号的拟合关系，解决了寻找标准弹道的问题。
3) 对标准弹道数据进行了进一步挖掘，求解出飞行时间的解析式，减少了迭代插值运算的次数，并且在制导过程中不再使用标准弹道文件，增加了对标准弹道的线下利用率，减少线上计算量，有效提升了算法效率。
4) 以处于滑翔段的高超声速飞行器为目标，改变射程、改变航路捷径、加入气动系数不确定性进行多次拦截仿真，验证了拦截弹拦截能力及预测制导方法的有效性。

参考文献

[1]	ZARCHAN P. Tactical and strategic missile guidance[M]. 6th ed. Reston: AIAA, 2012.
[2]	KUMAR P, DWIVEDI P N, BHATTACHARYA A, et al. Variable gain predictive PN guidance for interception of high speed re-entry targets[J]. IFAC-PapersOnLine, 2018, 51(1): 64-71. DOI:10.1016/j.ifacol.2018.05.012
[3]	李辕, 赵继广, 白国玉, 等. 基于预测碰撞点的剩余飞行时间估计方法[J]. 北京航空航天大学学报, 2016, 42(8): 1667-1674. LI Y, ZHAO J G, BAI G Y, et al. Method of time-to-go estimation based on predicted crack point[J]. Journal of Beijing University of Aeronautics and Astronautics, 2016, 42(8): 1667-1674. (in Chinese)
[4]	李辕, 闫梁, 赵继广, 等. 顺轨拦截模式剩余飞行时间估计方法[J]. 航空学报, 2015, 36(9): 3082-3091. LI Y, YAN L, ZHAO J G, et al. The method of time-to-go estimation for head-pursuit engagement[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(9): 3082-3091. (in Chinese)
[5]	DHANANJAY N, GHOSE D. Accurate time-to-go estimation for proportional navigation guidance[J]. Journal of Guidance, Control, and Dynamics, 2014, 37(4): 1378-1383. DOI:10.2514/1.G000082
[6]	GHOSH S, GHOSE D, RAHA S. Unified time-to-go algorithms for proportional navigation class of guidance[J]. Journal of Guidance, Control, and Dynamics, 2016, 39(6): 1-18.
[7]	GRUBIN C, SHULTZ P, SOUFL R.A boost guidance scheme for following a trajectory profile and satisfying injection constraints[C]//Astrodynamics Guidance and Control Conference, 1964: 639.
[8]	SALAMA M, HANCOCK G, SHOUMAN A.The predicted interception guidance method[C]//25th AIAA Aerospace Sciences Meeting.Reston: AIAA, 1987: 126.
[9]	张华伟, 董茜, 王文灿, 等. 基于预测命中点的反弹道导弹拦截方法研究[J]. 弹箭与制导学报, 2007, 27(2): 196-199. ZHANG H W, DONG Q, WANG W C, et al. Rearch way of intercepting ballistic missile based on the forecasting hitting position[J]. Journal of Projectiles Rockets, Missiles and Guidance, 2007, 27(2): 196-199. DOI:10.3969/j.issn.1673-9728.2007.02.062 (in Chinese)
[10]	ZHANG X, LEI H M, LI J, et al.Ballistic missile trajectory prediction and the solution algorithms for impact point prediction[C]//Proceedings of 2014 IEEE Chinese Guidance, Navigation and Control Conference.Piscataway: IEEE Press, 2014: 879-883.
[11]	ZHAO Y T, HU Y A, ZHANG Y A.Design and simulation of an optimal guidance law for ship-air missile based on interception point prediction[C]//2011 Chinese Control and Decision Conference (CCDC).Piscataway: IEEE Press, 2011: 3647-3652.
[12]	HAHN P V, FREDERICK R A, SLEGERS N. Predictive guidance of a projectile for hit-to-kill interception[J]. IEEE Transactions on Control Systems Technology, 2009, 17(4): 745-755. DOI:10.1109/TCST.2008.2004440
[13]	DWIVEDI P N, BHALE P G, BHATTACHARYYA A, et al. Generalized state estimation and model predictive guidance for spiraling and ballistic targets[J]. Journal of Guidance, Control, and Dynamics, 2014, 37(1): 243-264. DOI:10.2514/1.60075
[14]	SONG E J, TAHK M J. Real-time midcourse guidance with intercept point prediction[J]. Control Engineering Practice, 1998, 6(8): 957-967. DOI:10.1016/S0967-0661(98)00041-0
[15]	ZHANG J, YOU L Q, CHEN W C.Boost-phase guidance with neural network for interception of ballistic missile[C]//International Conference on Control, Automation and Information Sciences.Piscataway: IEEE Press, 2015: 426-431.
[16]	陈万春. 一种基于多项式拟合法的拦截弹预测制导方法: CN107766967A[P].2018-03-06.CHEN W C.A predictive guidance law of interceptor based on polynomial fitting method: CN107766967A[P].2018-03-06(in Chinese).
[17]	尤刘球. 动能拦截弹目标信息估计与预测制导方法研究[D]. 北京: 北京航空航天大学, 2017: 1-2. YOU L Q.Research on target information estimation and predictive guidance[D].Beijing: Beihang University, 2017: 1-2(in Chinese).
[18]	韩英宏, 雷延花, 梁卓, 等. 带侧窗动能杀伤器直接力姿态控制[J]. 航天控制, 2015, 33(4): 51-55. HAN Y H, LEI Y H, LIANG Z, et al. Direct force attitude control of kinetic kill vehicle with side window[J]. Aerospace Control, 2015, 33(4): 51-55. DOI:10.3969/j.issn.1006-3242.2015.04.009 (in Chinese)
[19]	谢经纬, 陈万春. 大气层外拦截弹建模与攻防效能分析[J]. 北京航空航天大学学报, 2018, 44(9): 1826-1838. XIE J W, CHEN W C. Modeling of exo-atmospheric interceptor and effectiveness analysis for penetration and defense[J]. Journal of Beijing University of Aeronautics and Astronautics, 2018, 44(9): 1826-1838. (in Chinese)
[20]	盛永智, 陈万春, 孟璇. 反导预测拦截及防御问题研究[C]//2005中国飞行力学学术年会, 2005: 110-115. SHENG Y Z, CHEN W C, MENG X.Research on predictive interception and defense of antimissile[C]//2005 Proceedings of China Annual Conference on Flight Mechanics, 2005: 110-115(in Chinese).
[21]	李静琳, 陈万春, 闵昌万. 高超末段机动突防/精确打击弹道建模与优化[J]. 北京航空航天大学学报, 2018, 44(3): 555-567. LI J L, CHEN W C, MIN C W. Terminal hypersonic trajectory modeling and optimization for maneuvering penetration and precision strike[J]. Journal of Beijing University of Aeronautics and Astronautics, 2018, 44(3): 555-567. (in Chinese)
[22]	周慧钟, 李忠应, 王瑾枚. 有翼导弹飞行动力学[M]. 北京: 北京航空航天大学出版社, 1983: 46-47. ZHOU H Z, LI Z Y, WANG J M. Flight dynamics of winged missiles[M]. Beijing: Beihang University Press, 1983: 46-47. (in Chinese)
[23]	陈桂秀. 用程序求解最小二乘拟合多项式的系数[J]. 青海师范大学学报(自然科学版), 2010, 26(3): 14-17. CHEN G X. Solve the least square curve fitting polynomial coefficient with program[J]. Journal of Qinghai Normal University(Natural Science Edition), 2010, 26(3): 14-17. DOI:10.3969/j.issn.1001-7542.2010.03.004 (in Chinese)
[24]	刘瑶, 张占月, 黄梓辰, 等. 基于拦截纵深的中段反导武器部署研究[J]. 指挥与控制学报, 2017, 3(2): 119-126. LIU Y, ZHANG Z Y, HUANG Z C, et al. Deployment of midcourse anti-missile weapon based on intercepting depth[J]. Journal of Command and Control, 2017, 3(2): 119-126. DOI:10.3969/j.issn.2096-0204.2017.02.0119 (in Chinese)
[25]	刘芳, 陈万春. PAC-3拦截弹六自由度反导建模与拦截仿真分析[J]. 飞行力学, 2012, 30(5): 440-443. LIU F, CHEN W C. PAC-3 interceptor 6DOF antimissile modeling and intercept simulation analysis[J]. Flight Dynamics, 2012, 30(5): 440-443. (in Chinese)