掩蔽库是飞机的重要防护工事, 目前通常使用侵彻类战斗部打击掩蔽库飞机^[1], 但随着配筋超高强混凝土、钢纤维增强混凝土、分层结构等技术应用^[2-4], 掩蔽库抗冲击侵彻能力有了较大提升, 侵彻类战斗部打击掩蔽库飞机的难度急剧增加。因此, 发展打击掩蔽库飞机的新手段变得日益重要。弹道式杀爆弹具有杀伤范围广、射程远、突防能力强、飞行时间短等优点, 如能用于毁伤掩蔽库飞机对提升对敌方机场飞机威慑能力有重要意义。杀爆弹打击掩蔽库飞机尚缺少相关理论与试验研究, 因此, 本文通过仿真评估的方法讨论这种新型打击方式的经济性与可行性。

对于杀爆弹打击飞机目标的试验与仿真研究, Moxnes等仿真分析了杀爆弹破片的飞散特性^[5], 陈新祥提出使用射击迹线法评估破片对空中目标的毁伤效果^[6], 王国庆等提供了一种破片式战斗部打击飞机目标的三维建模与毁伤评估方法^[7]。然而, 掩蔽库对破片遮挡严重, 飞机中弹数少, 随机涨落大, 目前方法均需在建模精度与弹目交会计算速度之间进行取舍, 难以满足毁伤评估所需的大规模蒙特卡罗仿真精确的需求。

与传统火力筹划中耗弹量与武器系统精度为单变量映射关系不同^[8], 杀爆弹打击掩蔽库飞机时, 最小耗弹量与起爆设置、终点弹道参数、瞄准点位置等多变量存在高度的非线性耦合, 而且与传统单瞄准点打击方案不同^[7], 由于掩蔽库遮挡, 单个瞄准点打击掩蔽库飞机易出现顾此失彼、毁伤不均等现象。

针对上述问题, 本文以某典型飞机与单机、双机掩蔽库为研究对象, 通过一/二次曲面拟合的方法构建目标三维曲面模型, 并建立破片的三维飞散模型, 结合弹目坐标系转换关系提出了弹目交会算法, 基于建立的易损性模型提出了飞机毁伤评估算法; 然后, 提出了基于网格法的飞机平均毁伤概率计算方法, 以及以整体毁伤效果最优为目标的瞄准点位置与耗弹量寻优方法, 并开展了打击方案优化方法研究; 最后, 仿真评估了某弹道式杀爆弹在不同落速、落角、爆高、交会角、CEP工况下重度毁伤掩蔽库飞机所需的最小耗弹量, 给出了打击方案优化策略。

1 飞机与掩蔽库的三维建模

以图 1的飞机为研究对象，该机机长15.3 m，翼展10.9 m，高度5.3 m^[9]。为兼顾建模精度与计算速度, 采用一/二次曲面拟合方法建立由圆锥、圆柱、三角形等10个曲面构成^[10]的该机简化模型。

以图 2飞机掩蔽库为对象, 利用一/二次曲面拟合减少曲面数量的思路, 将单机掩蔽库简化为由圆柱面、环形面等3个曲面构成。对单机掩蔽库结构进行合理拓展得到由11个曲面构成的双机掩蔽库。

假设飞机存放在掩蔽库正中间, 建立飞机系下战斗机存放模型见图 3~4。

2 杀爆战斗部破片的三维飞散模型

杀爆战斗部可简化为若干个圆台面, 建立爆点直角坐标(b), 见图 5。破片位置可依据战斗部圆台尺寸进行紧密排布得到。

弹道导弹在终点弹道的攻角与侧滑角一般近似为零, 根据图 5, 考虑弹头落速叠加, 得到任意破片e在爆点直角坐标系下的初速为^[10]

(1)

式中: v_s表示破片的抛撒速度; γ为破片的方位角在[0, 2π]区间服从均匀分布; ϕ为破片的抛射角, 在[ϕ_min, ϕ_max]区间服从截短的正态分布^[11]; v_m表示弹头落速。

破片轨迹近似为直线, 速度衰减规律为^[12]

(2)

式中: v₀表示破片初速; s表示飞行距离; k表示速度衰减系数; v表示破片飞行距离s后的速度。

考虑到掩蔽库对冲击波遮挡严重、瞄准点散布范围较大、爆高较高等因素, 本文不计算作用距离较短的冲击波超压对飞机的毁伤。

3 弹目交会计算 3.1 主要坐标系的建立与相互转换关系

主要坐标系见图 6, 定义如下:

1) 爆点直角坐标系(b): 原点位于弹尖, z轴沿弹轴向后, y轴位于射面向上, x轴由右手定则确定;

2) 瞄准坐标系(Aim): 原点位于弹头落速与地面的交点, x轴平行于弹头落速在地面的投影, z轴指向天, y轴由右手定则确定;

3) 飞机坐标系(A): 原点位于机头, y轴沿机身轴线指向后, z轴指向天, x轴由右手定则确定。

图 6中: 爆高H为弹尖距地面高度; 落角α为弹头落速与地面夹角; 交会角β为x_Aim轴与-y_A轴夹角。

根据图 6, 爆点直角坐标系至瞄准坐标系的坐标转换矩阵为

(3)

根据图 6, 瞄准坐标系至飞机坐标系的坐标转换矩阵为

(4)

根据图 6, 起爆点O_b在瞄准系下坐标为

(5)

根据图 6, 瞄准坐标系原点O_Aim在飞机坐标系下的坐标为

(6)

式中, x_Hit, y_Hit为瞄准点在飞机系下x, y坐标。

进而得到, 对于任意破片e, 其坐标由爆点直角坐标系转换至飞机坐标系的方法为

(7)

对于任意破片e, 其速度由爆点直角坐标系转换至飞机坐标系的方法为

(8)

3.2 弹目交会计算

弹目交会计算流程如下: 首先, 将所有破片位置与速度由爆点直角坐标系转换至飞机坐标系, 逐一计算各破片与所有曲面的交会情况^[12]; 然后, 假定破片只能与至多一个曲面交会, 基于最小飞行距离原则确定破片与飞机各曲面、掩蔽库各曲面以及地面的最终交会情况, 进而得到与飞机交会破片的着靶位置、速度、角度等参数。

3.3 飞机的毁伤计算

图 1中飞机要害舱段的位置分布见图 7。图 7中, 油箱、发动机、驾驶舱的破坏可引起飞机Ⅰ级毁伤(重度毁伤)^[13], 建立Ⅰ级毁伤树见图 8。

破片b对要害舱段i的击穿概率为^[13]

(9)

式中: E_T=4.5×10⁸ J/m³为撞击比动能阈值; e_b为破片b撞击比动能(J/m³), 可表示为

(10)

式中: m_b, v_b, θ_b, Φ_b分别为破片b的质量(kg)、着靶速度(m/s)、着靶角度、形状系数(m²/kg^2/3); h_d为等效靶厚度(m)。

破片毁伤包括击穿、引燃、引爆等, 但破片引燃油箱、引爆弹药概率极小, 因此仅考虑击穿毁伤效应, 即, 破片对要害舱段i毁伤概率为^[10]

(11)

式中: N_i为与要害舱段i撞击的破片数量; r_i为要害舱段i的易损性系数。

由于缺乏图 7飞机要害舱段材料、结构数据, 其等效靶厚度、易损性系数参照文献[13-14]。

根据(11)式易证, M发杀爆弹联合打击条件下, 要害舱段i的毁伤概率为

(12)

式中, P_D^{i, k}为杀爆弹k对要害舱段i的毁伤概率。

根据图 8毁伤树, 得到飞机毁伤概率为

(13)

式中, P_DE, P_DC, P_DB, P_DW分别为发动机、驾驶舱、机身整体油箱、机翼油箱的毁伤概率, 计算方法见(11)~(12)式。

据此得到图 1飞机毁伤概率计算流程如下:

利用弹目交会模型计算得到与各要害舱段撞击破片群的数量、着靶速度/角度/位置, 基于对图 7分析得到的各要害舱段易损性系数, 结合(9)~(11)式计算破片群各要害舱段的击穿毁伤概率, 再根据(13)式即可得到飞机毁伤概率。

若M发杀爆弹起爆高度、终点弹道参数、瞄准点位置相同, 且单发弹毁伤概率为P_RS, 则联合打击条件下, 飞机的毁伤概率为

(14)

若飞机重度毁伤概率阈值为P_DT, 则起爆高度、终点弹道参数、瞄准点位置相同时, 杀爆弹最小耗弹量为

(15)

式中, ⌈⋅⌉表示对实数“·”向上取整。

4 基于最小耗弹量的打击方案优化

经济性与可行性是杀爆弹能否用于打击掩蔽库飞机的核心问题, 即是否存在耗弹量可接受的打击方案, 因此就如下2个优化问题开展研究:

1) 终点弹道参数确定条件下, 杀爆弹最小耗弹量的寻优计算。弹道导弹射程远、末端机动少, 落点位置小幅变化对落速、落角、交会角、CEP等终点弹道参数影响较小; 另一方面, 打击固定地面目标时, 杀爆弹通常采用固定高度起爆^[15], 因此, 当弹道方案、起爆设置、导航体制确定时, 导弹在目标附近任意位置的终点弹道参数可近似认为是定值, 即, 终点弹道参数由弹道方案、起爆设置、导航体制的选择决定。而根据弹目交会模型可知, 终点弹道参数确定时, 杀爆弹的瞄准点位置决定了目标的毁伤效果。因此, 杀爆弹最小耗弹量计算可等效为如下问题: 以重度毁伤所有飞机为目标, 耗弹量与瞄准点位置的寻优。

2) 打击方案优化。基于终点弹道参数与最小耗弹量的映射关系, 以最小耗弹量为指标, 开展终点弹道参数的寻优分析, 考虑弹道方案、起爆设置、导航体制等约束, 提出打击方案优化策略, 并评估经济性与可行性。

为了简化问题, 本文以多发杀爆弹在相同弹道方案、起爆设置、导航体制下对目标进行打击为研究背景, 讨论打击方案优化方法。

4.1 终点弹道参数与最小耗弹量映射关系的计算方法 4.1.1 瞄准点对应的飞机平均毁伤概率计算

瞄准点位置寻优需要进行大量位置搜索, 受落点及破片随机散布影响, 如果各位置均开展大量的弹目交会蒙特卡罗计算, 计算耗时极长。因此, 本文采用对目标区域网格化, 将各网格点作为备选瞄准点的思路以大幅减少弹目交会计算。考虑到数量众多的网格点对应的飞机平均毁伤概率计算量仍极大, 而且CEP发生变化时, 经典方法^[16]需对每个网格点根据落点散布重新计算飞机平均毁伤概率。对此, 本文提出先不考虑CEP影响计算各网格点对应的飞机毁伤概率, 然后考虑CEP影响利用各网格点平均毁伤概率的加权平均值代替蒙特卡罗仿真的毁伤概率平均值, 从而避免CEP变化时各网格点对应飞机平均毁伤概率的重新计算。基于以上分析, 不考虑CEP影响时各网格点对应飞机平均毁伤概率计算方法如下:

step1  初始化落速、落角、爆高、交会角、CEP等终点弹道参数, 建立以掩蔽库为中心, 边长为L_SM正方形区域, 生成N_M×N_M个网格点;

step2  遍历x方向网格序号j_x=1→N_M;

step3  遍历y方向网格序号i_y=1→N_M;

step4  计算网格点(j_x, i_y)对应飞机系坐标

(16)

式中: x_center, y_center分别为掩蔽库中心在飞机坐标系的x, y坐标, 经核算, x_center=0 m, y_center=8 m。

step5  考虑到破片飞散随机性的影响, 杀爆弹以网格点(j_x, i_y)为瞄准点进行N_S次(N_S表示蒙特卡罗仿真次数)弹目交会计算, 统计各次仿真中飞机1、飞机2各要害舱段的毁伤概率, 循环结束后转入step6;

step6  计算网格点(j_x, i_y)对应飞机1、飞机2要害舱段平均毁伤概率(以飞机1发动机为例)

(17)

式中: P_{E, 1}(k)为飞机1在第k次仿真中发动机的毁伤概率。

通过以上步骤得到不考虑CEP影响时各网格点对应飞机各要害舱段平均毁伤概率后, 考虑CEP影响时网格点(j_x, i_y)对应飞机1的各要害舱段平均毁伤概率计算方法如下(以发动机为例)

(18)

式中, ，Δ为CEP取值。

杀爆弹以网格点(j_x, i_y)为瞄准点打击掩蔽库飞机1的平均毁伤概率可依据(13)式计算得到。

飞机2的平均毁伤概率也按上述步骤计算。

4.1.2 最小耗弹量寻优方法

以全部飞机整体毁伤效果最优为目标, 提出如下适应值函数

(19)

式中, P_R表示各飞机平均毁伤概率组成的向量。

显然, 终点弹道参数确定时, 适应值由瞄准点位置与耗弹量决定, 且适应值越大, 整体毁伤效果越好。若适应值函数满足

(20)

则所有飞机全部重度毁伤, 其对应的瞄准点位置与耗弹量为可行解。在适应值函数满足(20)式前提下, 使耗弹量最小, 则其对应的瞄准点位置与耗弹量为最优解。

1) 瞄准点只有1个的情况

根据(15)式, 单瞄准点条件下杀爆弹联合打击单机掩蔽库飞机时, 最小耗弹量与最优瞄准点位置计算方法为

(21)

根据(15)式, 单瞄准点条件下杀爆弹联合打击双机掩蔽库飞机时, 最小耗弹量与最优瞄准点位置计算方法为

(22)

(21)~(22)式可通过穷举法寻优, 算法只需遍历M=N_M×N_M个网格点, 计算量尚可接受。

显然, N_SA是耗弹量一个保守的上限。如果单瞄准点条件下的最小耗弹量大于1个, 则需开展多瞄准点条件下的最小耗弹量寻优。

2) 瞄准点数量大于1个的情况

多瞄准点条件下, 需要对每个导弹的瞄准点进行独立寻优, 当导弹数量为n, 穷举法的计算时间复杂度为O(Mⁿ), 在网格数量大, 导弹数量多的情况下, 穷举法计算量将十分巨大。

对此, 本文利用改进粒子群算法对瞄准点位置与耗弹量进行寻优, 步骤如下:

step1  遍历导弹个数n=2→N_SA;

step2  生成L个粒子的初始位置与速度^[17]

(23)

式中: x_i^C, v_i^C表示粒子i的位置与速度, 均为2n×1维向量, 且v_i^C服从(0, L_SM/10)正态分布, x_i^C由n发杀爆弹瞄准点的x, y坐标构成

(24)

x_Hit(k), y_Hit(k), v_x^C(k), v_y^C(k)分别为第k个杀爆弹瞄准点在飞机系A下的x, y坐标及对应粒子速度, 且x_Hit(k)服从区间[x_center-L_SM/2, x_center+L_SM/2]的均匀分布, y_Hit(k)服从区间[y_center-L_SM/2, y_center+L_SM/2]的均匀分布。

为加快算法收敛速度, 本文提出如下措施:

1) 第1个粒子由上一轮优化得到的n-1个最优瞄准点坐标与单个瞄准点的最优坐标构成;

2) 占总数1/5的粒子位置的前n-1个瞄准点坐标由上一轮优化得到的最优瞄准点坐标构成;

step3  将粒子个体最优位置初始化为当前位置, 即P^C(0)=X^C(0), 初始化各粒子适应值为0。

step4  循环代数t=1→T_max(最大代数);

step5  计算各粒子适应值(以粒子i为例)。

首先, 对于粒子i的第k发杀爆弹, 根据最近距离法确定其对应的网格点序号(j_x, i_y)为

(25)

如序号(j_x, i_y)超出1→N_M, 即粒子位置超边界, 则杀爆弹k对所有飞机平均毁伤概率均为0。

然后, 参照(18)式、依据(13)式计算粒子i各杀爆弹对各飞机的各要害舱段、飞机整体的平均毁伤概率;

最后, 根据(12)式、(13)式计算粒子i所有杀爆弹联合打击掩蔽库各飞机的毁伤概率, 再根据(19)式得到粒子i的适应值。

step6  更新各粒子的个体与全局最优位置

(26)

(27)

式中: F_fit(x_i^C(t))为粒子位置x_i^C(t)处适应值。

如果最优位置的适应值满足(20)式, 则杀爆弹能重度够毁伤全部飞机, 输出各导弹最优瞄准点位置p_g^C(t)与最小耗弹量n, 并结束计算。

step7  更新各粒子的速度与位置

(28)

式中: c₁, c₂为学习因子, 非负常数, 默认为2;r₁, r₂是在[0, 1]内均匀分布的随机数。

定义速度限制函数

(29)

取v_min^C=L_SM/2, v_max^C=-L_SM/2, 更新粒子速度

(30)

更新各粒子的位置

(31)

step8  超边界粒子的反弹处理。粒子位置若处于网格边界外, 飞机的毁伤概率便为0, 由于惯性机制, 粒子将长时间处于网格边界外, 这将导致大量无效寻优。对此, 本文提出了粒子反弹策略: 粒子位置超越边界时, 粒子速度立即反向。以瞄准点x坐标为例, 说明反弹处理策略

(32)

显然, 粒子反弹策略确保了粒子位置始终处于网格边界范围内, 提高了搜索效率。

4.2 打击方案优化方法

弹道导弹的弹道方案、起爆设置、导航体制决定了落速、落角、爆高、交会角、CEP等终点弹道参数, 因此, 打击方案优化步骤如下:

step1  根据弹道方案、起爆设置、导航体制的调整范围得到终点弹道参数取值集合;

step2  基于终点弹道参数与最小耗弹量映射关系, 分析终点弹道参数对最小耗弹量的影响;

step3  比较不同弹道方案、起爆设置、导航体制下的最小耗弹量, 提出打击方案优化策略。

5 仿真案例分析

杀爆弹参数: 破片数量为10 000枚, 破片初速为2 000 m/s, 直径10 mm, 密度17.6×10³ kg/m³, 速度衰减系数为1×10^-3m^-1, 抛射角为30°~90°。其他仿真参数: 重度毁伤阈值P_DT=0.95, 蒙特卡罗仿真次数N_S=100, 计算网格边长L_SM=100 m, 各方向网格点个数N_M=101;粒子群算法的粒子数L=100, 最大代数T_max=10 000。

仿真实例如下: 以落速600 m/s、落角70°、爆高20 m、交会角30°、CEP为10 m为基准案例, 计算弹道式杀爆弹毁伤单机掩蔽库、双机掩蔽库中所有飞机的最小耗弹量, 工况如下:

工况1, 仅落速在500~1 200 m/s内变化; 工况2, 仅落角在30°~90°内变化; 工况3, 仅爆高在10~40 m内变化; 工况4, 仅交会角在0°~90°内变化; 工况5, 仅CEP在5~20 m内变化。

针对基准案例仿真得到杀爆弹打击单机、双机掩蔽库飞机毁伤过程见图 9~10，图中红色、绿色破片分别为已撞击、飞行破片。可见, 一/二次曲面拟合方法不仅确保了建模精度, 而且通过大幅减少曲面数量提高了弹目交会计算速度, 有利于开展大规模蒙特卡罗仿真^[18]。

对基准案例仿真得到杀爆弹打击双机掩蔽库飞机1毁伤概率空间分布特性见图 11~12。可见, 毁伤概率比较大的区域在y≥18和y≤-11, 与图 4掩蔽库对破片遮挡效果相符; CEP对飞机毁伤概率分布有较强的展宽与压低效果, 体现了平均效应的影响, 验证了网格法计算飞机毁伤概率的准确性。另一方面, 与经典方法^[16]对每个网格点单独进行落点随机散布分析计算目标平均毁伤概率相比, 本文算法综合利用所有网格数据计算目标平均毁伤概率不仅数据利用充分、计算准确, 而且大幅减少了弹目交会计算, 提升了计算效率。

根据工况1~5, 仿真得到最小耗弹量随各终点弹道参数变化情况见图 13~17。

根据图 13~17, 可得到如下结论:

1) 落速越大, 最小耗弹量越少, 这是由于增加落速可以增大破片动能、降低前向盲区;

2) 随着落角的增加, 最小耗弹量呈现先减后增的变化, 这是由于落角过大时, 掩蔽库对破片遮挡严重, 破片对飞机的覆盖面积很小, 而落角过小时, 起爆点需要距离飞机很远才能实现破片对飞机的有效覆盖, 导致单位面积的中弹数很少;

3) 爆高越小, 最小耗弹量越少, 这是由于降低爆高可以有效增加单位面积的破片密度而且可以增加飞机的迎弹角度范围, 提高飞机的中弹数, 从冲击波超压毁伤角度来看, 低爆高还有助于提高作用距离较短的冲击波超压毁伤飞机的概率;

4) 交会角越小, 最小耗弹量越少, 这是由于减少交会角可以增加飞机的迎弹面积, 减少掩蔽库对破片的遮挡;

5) CEP越小, 最小耗弹量越少, 这是由于减小CEP可降低小概率毁伤区对飞机中弹数的影响;

6) 爆高与CEP对最小耗弹量的影响最为显著, 因此采用弹道式杀爆弹打击掩蔽库飞机的策略首先是降低爆高、提高打击精度, 其次是在发射阵地可选时, 射向与掩蔽库轴线的夹角(交会角)应尽可能小, 在飞行弹道可调时, 根据杀爆弹破片飞散特点优化弹道倾角, 并酌情增加落速;

7) 使用组合打击方案的双机掩蔽库最小耗弹量与单机掩蔽库的最小耗弹量相当, 可见, 组合打击方案可有效降低打击双机掩蔽库飞机成本。

6 结论

本文针对弹道式杀爆弹打击掩蔽库飞行开展了打击方案优化研究, 并得到如下结论:

1) 提出了基于一/二次曲面拟合的飞机、单机掩蔽库、双机掩蔽库的三维建模方法, 该建模方法在确保建模精度的前提下, 有效减少了建模曲面数量, 可实现弹目交会的快速计算, 有利于开展大规模蒙特卡罗仿真研究;

2) 给出了杀爆战斗部破片速度与位置的三维建模方法, 构建了杀爆弹打击掩蔽库飞机的坐标转换关系模型, 提出了弹目交会算法;

3) 提出了基于网格法的飞机平均毁伤概率计算方法, 减少了CEP变化时的计算量;

4) 提出以全部飞机整体毁伤效果最优为目标的多弹联合打击掩蔽库飞机的最小耗弹量与最优瞄准点位置寻优方法;

5) 改进了粒子群算法, 提出了粒子反弹策略, 即粒子位置超出计算边界时粒子速度立即反向, 避免了粒子进入无效搜索区, 提高了寻优效率;

6) 开展了弹道式杀爆弹打击掩蔽库飞机的仿真研究, 仿真结果表明, 爆高与CEP对最小耗弹量的影响最为显著, 因此, 弹道式杀爆弹打击掩蔽库飞机的策略重点是降低爆高、提高打击精度, 其次是减小交会角, 优化落角, 增加落速;

7) 采用多瞄准点组合打击方案可以有效降低杀爆弹打击掩蔽库飞机的耗弹量;

8) 多瞄准点组合打击双机掩蔽库飞机的寻优方法可拓展至弹道式杀爆弹对多个掩蔽库、多架飞机的组合打击方案优化。