多约束条件下多目标反导最优分配策略研究

马长李，张 涛，李 振

（海军研究院，北京，100161）

0 引 言

弹道导弹中段拦截是指在大气层外将来袭弹道导弹摧毁于自由段［1］。中段拦截具有弹道预报难度低、响应时间长、直接力控制简便等优点，是最为安全、可靠的弹道导弹防御方式。

对于中段拦截问题，最为重要的便是射前诸元解算，包括发射时刻、拦截时刻和发射方位角等。上述诸元直接影响了拦截弹的弹道形态以及末端交战场景，对拦截效能影响深远。同时，考虑拦截弹自身作战指标（如射程和射高），发射窗口和拦截弧段也会对拦截效能产生影响。高度动态变化的战场态势还对诸元解算方法的快速性提出了较高的要求［2］。因此，发射诸元解算是一个集精度、约束和速度于一体的高度非线性问题。文献［2］采用最小二乘迭代的方式计算发射诸元，但其主动段飞行模式固定，不能适应目标机动的情况，且计算耗时很长；文献［3］在二维平面内对拦截时间窗口问题进行了研究，但模型过于简单，缺乏实际意义；文献［4］考虑地基雷达预警时间，建立了反导预警时间模型，但计算中采用的是目标平均速度，计算精度不高。

求解拦截弹发射诸元是制定拦截策略的前提，对于在实际应用中经常出现的多拦截器-多目标的场景，还必须针对各目标的飞行状态和拦截弹的部署情况（如位置、数量和作战能力等），将各目标分配给各拦截弹。上述分配的依据便是发射诸元计算结果、拦截弹作战指标和目标威胁程度，而分配的过程则需要通过优化算法来完成。多目标拦截任务分配问题属于武器-目标分配范畴［5-7］，在作战决策和火力运用中有重要地位。国内外学者对武器系统、编队协同作战及舰艇编队防空等场景中的目标分配模型进行了大量研究［8-11］。

本文首先以二体理论和Gauss 问题为基础，以发射时刻和命中时刻为设计变量，提出了一种考虑多约束条件的拦截弹发射诸元快速解算方法，解决了一对一条件下的发射窗口、拦截弧段和发射诸元快速求解问题。其次，建立了多对多最优分配模型，包括性能指标、目标威胁指标和拦截弹作战指标等。另外，设计了一种改进的粒子群优化算法，结合上述优化模型对最优分配策略进行求解。通过引入动态种群规模和种群变异策略，提高了粒子群算法的搜索效率，保证了最优分配策略的快速生成。最后通过仿真计算验证了方法的准确性和可行性。

1 发射诸元快速解算方法

本文定义的发射诸元包括发射方位角、发射时刻和命中时刻，同时定义发射窗口为所有发射时刻的集合。每个发射时刻都对应目标弹道上的一个弧段，记整个发射窗口对应的所有弧段的集合为拦截弧段。

该方法能够根据目标实际运动状态快速求出拦截弧段，预测命中点可通过开普勒轨道递推方程快速求解，Gauss 问题通过几次牛顿迭代即可求解。因此，对于目标进行轨道机动情形，采用本文提出的方法可以针对机动情形快速求出发射窗口和拦截弧段。

首先作如下设定：

a）目标轨迹已由预警探测系统预报得出，为简化计算复杂度，这里暂不引入目标机动项。

b）设红方和蓝方分别为攻方和拦截方。

记仿真初始时刻T0=0，此时红方所有目标已被蓝方准确探测跟踪，蓝方具备诸元解算条件。

设蓝方完成解算、分发、发射等准备工作所需的时间为T1，记发射时刻为T2，显然T2≥T1。

c）记发射点在地心惯性系下的位置矢量为r0，经纬度分别为λ0和φ0。

本文设计发射诸元计算方法如下：

a）给定T2的初值，取T2=T1。

b）给定T3的初值。估计拦截弹从发射到命中的飞行时间为TM，取T3=T2+TM。

c）根据弹道预报结果，求出T3时刻目标在地心惯性系下的位置矢量rT（即预测命中点）、速度矢量vT、经度λT和纬度φT。

d）根据r0、rT和TM，通过求解Gauss 问题得到预测命中点和发射点间的椭圆轨道，同时得到拦截弹终端速度矢量vM和拦截交会角。

e）判断如下条件：

1）Gauss问题有解；

2）预测命中点高度不超过拦截弹发射高度范围；

3）预测命中点到拦截弹发射点的射程L0T不超过拦截弹射程范围。

若以上条件都满足，则输出T2和T3作为可行解，否则令T3=T3+ Δt3并返回步骤2，其中Δt3为命中时刻的搜索步长；若T3达到搜索上限（预测命中点高度低于拦截弹射高下限）仍没有满足条件的解，令T2=T2+ Δt2并返回步骤1，其中Δt2为发射时刻的搜索步长；若T2达到搜索上限后仍无可行解，说明该拦截弹不具备对该来袭导弹的拦截条件。

拦截交会角［12］为弹目速度矢量夹角的补角，其计算方式如下：

射程的计算方式如下：

式中Re为地球平均半径；μ0T的计算方式如下：

综上，记发射窗口和拦截弧段分别为Ωlaunch和Ωintercept，则

上述方法能够快速确定某发射点处的拦截弹对某目标的发射窗口和拦截弧段。在Ωlaunch中选择具体的发射时刻T2，在Ωintercept中选择具体的拦截时刻T3，则发射方位角的计算方式如下：

2 多对多最优分配模型

设有拦截弹m个，来袭目标n个；规定每个拦截弹只能拦截一个目标，每个目标可分配给任意一个拦截弹。多目标分配的优化模型为

式中Vij为设计决策变量，当第i个拦截弹不拦截第j个目标时取0，当第i个拦截弹拦截第j个目标时取1；Wj为第j个目标的威胁程度；Pij为第i个拦截弹对第j个目标的拦截效能。

2.1 目标威胁指标

目标威胁程度与目标特性、拦截难度和预计落点等密切相关。

2.1.1 目标特性

目标特性主要指目标红外特性、电磁特性和空间分布等信息，其关键定量因素主要是温度、运动状态、雷达截面积（Radar Cross Section，RCS）等。

a）温度。

在大气层外飞行时，真弹头和诱饵受光照影响会吸收热量并升温。首先给出比热容计算公式：

式中CQ为比热容；ΔQ为吸热量；m为飞行器质量；ΔT为物体温度变化。

设真弹头和诱饵材料的比热容相近，由于真弹头的质量一般远大于诱饵质量，根据式（7）可知在吸收相同热量的条件下，诱饵的温度变化率远大于真弹头，因此可将目标温度变化率作为判断诱饵与真弹头的依据。取温度特性函数为

其中，

式中x1为目标温度的变化量；b1，a1分别为温度变化的上下界。根据文献［5］，此处取b1=180 K，a1=230 K。

b）运动状态。

由于诱饵一般无姿轨控设备，受分离干扰力和干扰力矩的影响，其在与母舱分离前后的速度变化非常明显。而真弹头质量大且具有姿轨控设备，其在与母舱分离前后的速度变化相对较小。

考虑目标运动状态探测过程的时间延迟、分离力大小和真弹头/诱饵的质量规模，本文取运动特性威胁度函数为

式中x2为目标运动特性的变化量；b2，a2分别为运动状态变化判断的上下界。对于位置状态来说，取a2=30 m，b2=40 m；对于速度状态来说，取a2=20 m/s，b2=35 m/s。

c）雷达散射截面积（Radar Cross Section， RCS）。

弹道导弹的中段运动特征是真假弹头识别所依据的主要特征之一。真弹头有姿控系统，其RCS数据序列受进动影响而呈现周期性；而诱饵一般无姿控机构，其旋转和翻滚可视作随机运动。RCS数据序列规律性较弱，变化幅度更大。

根据文献［5］，本文取关于RCS的特性威胁度函数为u3= exp(-x3)，其中x3为目标RCS变化量。

2.1.2 预计落点

通过弹道预报可知目标的打击对象。根据被打击对象的重要程度可评估目标威胁程度。设被攻击区域有M个，各区域按照重要程度排序。设某目标j攻击的是第ℓ个目标，取其落点威胁函数为

式中εℓ为第ℓ个蓝方目标的相对重要度，εℓ∈(0，1)。

获得u1～u4后，对各威胁函数加权求和。若某威胁度函数在各目标间的差异不大，则为其分配较小权重；否则说明该特性对目标识别明显，应为其分配较大权重。目标威胁指标计算方式如下：

2.2 拦截弹作战指标

本文设定拦截弹作战指标由发射窗口、拦截弧段和平均拦截交会角决定。

根据上文中的定义，本文认为发射窗口和拦截弧段越长则拦截概率越高。对于交会角，假设拦截弹末制导采用比例导引，则指令加速度大小为

式中k为导引系数；|ṙ|为弹目相对速度；ω为弹目视线角速度。

迎面拦截时弹目相对速度（即|ṙ|）很大，此时需求加速度也很大；反之，顺向拦截时需求加速度相对较小。由交会角定义可知，交会角越大越利于末端拦截，因此认为大交会角条件下的拦截概率更高。

记第i个拦截弹对第j个目标的发射窗口宽度为|，拦截弧段宽度为，平均拦截交会角为；记所有|中的最大值为，所有中的最大值为，所有中的最大值为。综上，拦截弹作战指标的计算方式如下：

3 武器目标最优分配算法

当有m个拦截弹和n个目标时，共有L=m⋅n种拦截方式。当L较小时，可以采用枚举方法获得最佳拦截方式；当L很大时，枚举法所产生的计算量将非常大，过长的计算时间将影响指挥控制系统的实时性和响应能力。多目标最优分配实际上是从有限种组合中选择一种的过程，优化过程中可以不考虑约束条件。因此，可以将目标最优分配问题转化为无约束优化问题，并采用数值优化算法求解。

求解优化问题的方法主要有两类：确定性方法和随机性方法。确定性方法收敛速度快，但对初值较为敏感。近年来，一些随机性的方法例如遗传算法、粒子群优化算法等已成功应用到优化问题中。其中，粒子群优化由于易操作、收敛快而受到广泛关注。

基本粒子群优化算法（Particle Swarm Optimization，PSO）是一种基于迭代模式的优化算法，搜索空间中的每个粒子都对应一个由优化函数所决定的适应度，以评价该粒子当前位置的优劣。设种群的规模为N，搜索空间为D维，则粒子群中第i个粒子的位置和速度为［13］

记粒子i迄今找到的最好位置为pi=（pi1，pi2，…,piD），所有粒子迄今找到的最好位置为pg=（pg1，pg2， …,pgD）。对于粒子i，其第j维的位置和速度变化如下：

式中c1，c2为学习因子；r1，r2为0～1 间的随机数；w为惯性权重；tpso为当前迭代次数。

较大的惯性权重可以加强算法的全局搜索能力，反之可以加强局部搜索能力。为提高算法搜索性能，采用指数递减的惯性权重：

式中w0为起始权重；wf为终端权重；d1，d2为控制因子；tmax为最大的迭代次数。

由于算法的随机性，会出现粒子聚集在某个或几个特定位置的情况，这主要取决于问题特性以及适应度函数的选择。为定量描述粒子状态，引入群体适应度方差的定义。设fi为第i个粒子的适应度，favg为群体平均适应度，则群体适应度方差为

其中，

群体适应度方差反映的是粒子群中所有粒子的聚集程度，δpso越小，则聚集程度越大，导致将使群体失去多样性，陷入早熟收敛，故当δpso

式中ηpso为服从标准正态分布的随机量。

4 仿真分析

设共有6 枚红方导弹对蓝方6 个地面目标进行攻击，共有拦截弹5 枚，故至少有1 个目标被遗漏，须通过优化分配使蓝方损失最小。

取粒子数为200，最大迭代次数为20。在计算发射诸元时，取T1=60 s，TM=300 s，Δt2=10 s，Δt3=30 s。其他条件如表1所示。

表1 拦截弹参数和目标初始状态Tab.1 The initial states of the missiles and the targets

采用式（13）计算各拦截弹对各目标的作战指标，得到结果如表2所示。使用本文提出的最优拦截策略制定方法，得到仿真结果如表3 所示。可以看出，在设定的攻防对抗场景下，来袭目标4和来袭目标6的威胁最大，应当首先拦截。以拦截弹5和目标1为例，给出拦截弹拦截弧段如图1和图2所示。

图1 拦截弹5对目标1的拦截弧段Fig.1 The interception arc of missile-5 with respect to target-1

图2 拦截弹5对目标1的拦截弧段Fig.2 The interception arc of missile-5 with respect to target-1

最后，通过仿真验证改进PSO算法对搜索能力的提升。采用基本PSO方法求解上述最优分配问题，其中取惯性权重为常值0.8，不采用式（20）所示的粒子扰动策略。由于PSO算法具有随机性，对基本PSO和改进PSO 均进行30 次独立仿真，并取各自最好的仿真结果进行对比，结果如图3所示。其中，适应度越小代表优化函数值越优，即结果的性能指标越好。仿真结果表明，本文提出的改进PSO方法较基本PSO方法而言，能够在更少的进化代数内获得更优的性能指标，搜索能力更强。

图3 基本PSO和改进PSO计算结果对比Fig.3 Comparison between the basic PSO and the improved PSO

5 结束语

针对弹道导弹中段拦截中的多对多作战场景，提出了一种目标最优分配算法。首先，基于诸元快速解算法构建了拦截效能与威胁指标的快速求解方法，在此基础上建立了目标分配优化模型。其次，提出了一种改进粒子群优化算法，通过种群变异策略提高了算法收敛速度，获得了更优的求解结果。最后，通过仿真验证了优化模型的正确性与求解算法的快速性。