基于改进粒子群算法的烟幕干扰弹投放策略研究 - 橙子薄荷的博客

如果不出意外的话, 这应该是最后一次参加数学建模比赛了. 这次参赛的体验很不错, 作为代码手, 能够顺利地从头到尾写完所有的程序、再画图、写论文中的算法部分, 我是非常非常满意滴! 虽然在建模中仍有一些遗憾, 论文手的发挥也不尽如人意, 但是我们完成了论文, 做出了很好的答案, 也不奢求更好的结果啦~

1 题目

烟幕干扰弹主要通过化学燃烧或爆炸分散形成烟幕或气溶胶云团，在目标前方特定空域形成遮蔽，干扰敌方导弹，具有成本低、效费比高等优点。随着烟幕干扰技术的不断发展，现已有多种投放方式完成烟幕干扰弹的定点精确抛撒，即在抛撒前能精确控制烟幕干扰弹到达预定位置，通过时间引信时序控制起爆时间。

现考虑运用无人机完成烟幕干扰弹的投放策略问题。具有长续航能力的无人机挂载某型烟幕干扰弹在特定空域巡飞，受领任务后，无人机投放烟幕干扰弹在来袭武器和保护目标之间形成烟幕遮蔽。每架无人机投放两枚烟幕干扰弹至少间隔 $1$ s。烟幕干扰弹脱离无人机后，在重力作用下运动。烟幕干扰弹起爆后瞬时形成球状烟幕云团，由于采用特定技术，该烟幕云团以 $3$ m/s 的速度匀速下沉。据试验数据知，云团中心 $10$ m 范围内的烟幕浓度在起爆 $20$ s 内可为目标提供有效遮蔽。

来袭武器为空地导弹，该型导弹飞行速度 $300$ m/s。导弹的飞行方向直指一个为掩护某半径 $7$ m、高 $10$ m 的圆柱形固定目标而专门设置的假目标。以假目标为原点，水平面为 $xy$ 平面，真目标下底面的圆心为 $(0,200,0)$。警戒雷达发现来袭导弹时，$3$ 枚导弹 $\text{M1}$、$\text{M2}$、$\text{M3}$ 分别位于 $(20000,0,2000)$、$(19000,600,2100)$、$(18000,-600,1900)$；$5$ 架无人机的位置信息分别为 $\text{FY1}(17800,0,1800)$、$\text{FY2}(12000,1400,1400)$、$\text{FY3}(6000,-3000,700)$、$\text{FY4}(11000,2000,1800)$、$\text{FY5}(13000,-2000,1300)$。

在导弹来袭过程中，通过投放烟幕干扰弹尽量避免来袭导弹发现真目标。控制中心在警戒雷达发现目标时，立即向无人机指派任务。无人机受领任务后，可根据需要瞬时调整飞行方向，然后以 $70\text{∼}140$ m/s 的速度等高度匀速直线飞行。每架无人机的航向、速度可不相同，但一旦确定就不再调整。

为实现更为有效的烟幕干扰效果，需设计烟幕干扰弹的投放策略，主要包括无人机飞行方向、飞行速度、烟幕干扰弹投放点、烟幕干扰弹起爆点等。请建立数学模型，针对不同情形，分别设计烟幕干扰弹的投放策略，使得多枚烟幕干扰弹对真目标的有效遮蔽时间尽可能长。不同烟幕干扰弹的遮蔽可不连续。

问题 1: 利用无人机 FY1 投放 1 枚烟幕干扰弹实施对 M1 的干扰，若 FY1 以 120 m/s 的速度朝向假目标方向飞行，受领任务 1.5 s 后即投放 1 枚烟幕干扰弹，间隔 3.6 s 后起爆。请给出烟幕干扰弹对 M1 的有效遮蔽时长。
问题 2: 利用无人机 FY1 投放 1 枚烟幕干扰弹实施对 M1 的干扰，确定 FY1 的飞行方向、飞行速度、烟幕干扰弹投放点、烟幕干扰弹起爆点，使得遮蔽时间尽可能长。
问题 3: 利用无人机 FY1 投放 3 枚烟幕干扰弹，实施对 M1 的干扰。请给出烟幕干扰弹的投放策略，并将结果保存到文件 result1.xlsx 中（模板文件见附件）。
问题 4: 利用 FY1、FY2、FY3 等 3 架无人机，各投放 1 枚烟幕干扰弹，实施对 M1 的干扰。请给出烟幕干扰弹的投放策略，并将结果保存到文件 result2.xlsx 中（模板文件见附件）。
问题 5: 利用 5 架无人机，每架无人机至多投放 3 枚烟幕干扰弹，实施对 M1、M2、M3 等 3 枚来袭导弹的干扰。请给出烟幕干扰弹的投放策略，并将结果保存到文件 result3.xlsx 中（模板文件见附件）。

2 问题分析

问题 1 中无人机的飞行速度、飞行方向，烟幕干扰弹的投放、起爆位置都是确定的，难点是“有效遮蔽”的判定。考虑对导弹与烟幕云团的位置关系作分类讨论，当导弹处于烟幕云团内部时，可判定为有效遮蔽；当导弹处于烟幕云团外部时，可以采用阴影锥法进行判定。以导弹为顶点、母线与烟幕云团相切形成“阴影锥”。当真目标圆柱体的“最小外接球”完全被包围在阴影锥里面，可判定为有效遮蔽。

问题 2 是一个单无人机单干扰弹单目标问题，在问题 1 数学模型的基础上，可延用构建的阴影锥模型进行有效遮蔽判定。本题最优投放策略与无人机飞行速度、无人机飞行方向、投放烟幕干扰弹时刻以及烟幕干扰弹起爆时刻等四个决策变量有关，据此列出的目标函数难以使用非线性规划等传统方法直接求解。考虑到模型解空间复杂度较高，存在大量局部最优解，可以先使用粒子群算法求出问题的模拟最优解，最后以数值算法求出较精确的数值解。

问题 3 是一个单无人机多干扰弹单目标问题，无人机和导弹状态与前两问相同，但投放的烟幕干扰弹数量增为 3 枚，则要考虑两枚烟幕干扰弹投放时间间隔至少为 1s。另外, 本题需改进阴影锥法，除了单个球状烟幕云团的有效遮蔽，还需考虑两个球状烟幕云团部分重合后产生的有效遮蔽。考虑到高维、非线性、不可微的情况，可以使用粒子群优化算法解决本题。

问题 4 是一个多无人机单干扰弹单目标问题，投放的烟幕干扰弹数量与问题 3 相同，但变为由三架无人机各投一枚，则两枚烟幕干扰弹投放时间间隔无特殊要求。有效遮蔽判定沿用问题 3 中改进的阴影锥法。本题难点在于三架无人机具有独立的起始位置、速度与航向角，约束条件更复杂，影响决策的变量增多且相互间耦合度低，在使用一般粒子群算法时计算时间较长、容易陷入局部最优解。因此应该使用有初始点的粒子群方法，先讨论简化问题下的特解，再将特殊粗略解作为一个初始粒子进入粒子群算法, 以加快求解速度、提高求解质量。

问题 5 是一个多无人机多干扰弹多目标问题，需要先将无人机和导弹建立对应的干扰关系。为了简化模型，我们不妨假设一架无人机至多遮蔽一个导弹、一个导弹可被多个无人机遮蔽的。先根据无人机与导弹的距离进行分组，分组情况可能为“1:1:3”（“1”代表 1 个导弹由 1 架无人机遮蔽）或“1:2:2”。简化为“1”时可以退化为问题 3 的情形，简化为“2”（1 个导弹由 2 架无人机遮蔽）或“3”（1 个导弹由 3 架无人机遮蔽）时没有已有的解决策略。因此我们需要重点解决“2”或者“3”的简化情况，可以使用粒子群优化算法分析，此时分析思路与问题 4 相似，需要使用初始粒子改进算法。

3 建模前的准备

3.1 导弹的运动

设导弹 $Mi(i = 1,2,3)$ 被警戒雷达发现时的坐标为 $(M_{i,x}(0),M_{i,y}(0),M_{i,z}(0))$ 。导弹固定朝假目标即坐标原点方向运动，运动方向用向量坐标表示为 $(-M_{i,x}(0), - M_{i,y}(0), - M_{i,z}(0))$ 。导弹运动速度固定为 $300m / s$ ，速度向量为：

\[\begin{array}{l} \vec {v} _ {M _ {i}} = (v _ {i, x}, v _ {i, y}, v _ {i, z}) \\ = \left(\frac {- 3 0 0 M _ {i , x} (0)}{\sqrt {M _ {i , x} ^ {2} (0) + M _ {i , y} ^ {2} (0) + M _ {i , z} ^ {2} (0)}}, \right. \\ \frac {- 3 0 0 M _ {i , y} (0)}{\sqrt {M _ {i , x} ^ {2} (0) + M _ {i , y} ^ {2} (0) + M _ {i , z} ^ {2} (0)}}, \\ \left. \frac {- 3 0 0 M _ {i , z} (0)}{\sqrt {M _ {i , x} ^ {2} (0) + M _ {i , y} ^ {2} (0) + M _ {i , z} ^ {2} (0)}}\right) \tag {1} \\ \end{array}\]

导弹 $Mi$ 在 $t$ 时刻的坐标为：

\[\begin{array}{l} P _ {M _ {i}} (t) = \left(M _ {i, x} (t), M _ {i, y} (t), M _ {i, z} (t)\right) \\ = (M _ {i, x} (0) - v _ {i, x} t, M _ {i, y} (0) - v _ {i, y} t, M _ {i, z} (t) - v _ {i, z} t) \tag {2} \\ \end{array}\]

图1 导弹、无人机及目标位置示意图

3.2 无人机、烟幕干扰弹和烟幕云团的运动

(1) 无人机的运动

无人机接收到指派任务后瞬间调整方向, 然后以 $70 \sim 140 \mathrm{~m} / \mathrm{s}$ 的速度等高度匀速飞行。假设初始时刻第 $j(j = 1,2,3,4,5)$ 架无人机的位置为 $(FYj_{x}(0), FYj_{y}(0), FYj_{z}(0))$ , $t$ 时刻第 $j$ 架无人机的位置如图1所示, 为:

\[P _ {F Y j} (t) = \left(F Y j _ {x} (0) + v t \cos \phi , F Y j _ {y} (0) + v t \sin \phi , F Y j _ {z} (0)\right) \tag {3}\]

(2) 烟幕干扰弹的运动

在 $t_{u}$ 时刻烟幕干扰弹被投放, 烟幕干扰弹脱离无人机作平抛运动, 水平方向速度和无人机运动速度相同, 竖直方向在重力的作用下作自由落体运动。在 $t$ 时刻, 烟幕干扰弹的位置为:

\[P _ {f a l l} (t) = \left(F Y j _ {x} (0) + v t \cos \phi , F Y j _ {y} (0) + v t \sin \phi , F Y j _ {z} (0) - \frac {1}{2} g \left(t - t _ {u}\right) ^ {2}\right) \tag {4}\]

(3) 烟幕云团的运动

烟幕干扰弹在投放后经过 $\Delta t$ 时间起爆, 形成半径 $r = 10 m$ 的烟幕云团, 烟幕云团水平方向的速度为 0 , 竖直方向以 $3 m / s$ 的速度作匀速直线运动。 $t$ 时刻烟幕云团的球心 $C$ 坐标为:

\[\begin{array}{l} P _ {C} (t) = \left(C _ {x} (t), C _ {y} (t), C _ {z} (t)\right) \\ = \left(F Y j _ {x} (0) + v \left(t _ {u} + \Delta t\right) \cos \phi , \right. \\ F Y j _ {y} (0) + v \left(t _ {u} + \Delta t\right) \sin \phi , \\ F Y j _ {z} (0) - \frac {1}{2} g \Delta t ^ {2} - 3 \Delta t \Bigg) \tag {5} \\ \end{array}\]

3.3 有效遮蔽判定

(1) 导弹在烟幕云团内部

此时导弹与真目标上任意一点的连线都与烟幕云团相交，即从导弹的视角出发，真目标被烟幕云团完全遮挡，判定为有效遮蔽，即：

\[\sqrt {\left(M _ {i , x} (t) - C _ {x} (t)\right) ^ {2} + \left(M _ {i , y} (t) - C _ {y} (t)\right) ^ {2} + \left(M _ {i , z} (t) - C _ {z} (t)\right) ^ {2}} \leq r \tag {6}\]

(2) 导弹在烟幕云团外部

如图2a所示，此时从导弹出发且与烟幕云团相切所有射线构成一个圆锥，我们称之为阴影锥。当真目标完全位于阴影锥内时，说明从导弹的视角出发，真目标被烟幕云团完全遮挡，判定为有效遮挡。

图2a 真目标包围球示意图

图2b 阴影锥遮挡模型示意图

对于真目标是否完全位于阴影锥的判定，我们采用简化处理，如图2b所示，用以真目标的中心 $T(0,200,5)$ 为球心， $R = \sqrt{7^2 + (\frac{10}{2})^2} = \sqrt{74}$ 的长度为半径的球体 $T$ ，将真目标圆柱刚好完全包裹。此时从导弹出发且与球体 $T$ 相切所有射线构成一个新圆锥，设阴影锥的半顶角为：

\[\alpha = \arcsin \left(\frac {1 0}{\Vert \vec {M _ {i} C}\Vert} \right) \tag {7}\]

新圆锥的半顶角为

\[\beta = \arcsin \left(\frac {\sqrt {7 4}}{\Vert \vec {M _ {i} T} \Vert }\right) \tag {8}\]

两个圆锥的轴的夹角为

\[\theta = \arccos \left(\frac {\vec {M _ {i}} C \cdot \vec {M _ {i}} T}{\Vert \vec {M _ {i}} C \Vert \cdot \Vert \vec {M _ {i}} T \Vert}\right) \tag {9}\]

若满足导弹在烟幕云团外：

\[\sqrt {\left(M _ {i , x} (t) - C _ {x} (t)\right) ^ {2} + \left(M _ {i , y} (t) - C _ {y} (t)\right) ^ {2} + \left(M _ {i , z} (t) - C _ {z} (t)\right) ^ {2}} > r \tag {10}\]

烟幕云团在导弹与真目标之间：

\[\Vert \vec {M _ {i}} T \Vert \geq \Vert \vec {C T} \Vert \tag {11}\]

且

\[\beta + \theta \leq \alpha \tag {12}\]

则说明此时球体 $T$ 被阴影锥完全包含，同时满足真目标圆柱被阴影锥完全包含。

4 问题1的建模与求解

4.1 问题1模型求解

运用问题1的已知条件直接进行求解。

Step1 计算 $t$ 时刻各物体的运动轨迹函数

1）无人机FY1的运动轨迹

已知无人机FY1速度 $v_{FY1} = 120m / s$ ，初始位置 $P_{FY1}(0) = (17800,0,1800)$ ，目标位置 $P_{FY1,end} = (0,0,1800)$ ，则 $\phi = 180^{\circ}$ 。代入(3)得FY1的位置函数：

\[P _ {F Y 1} (t) = (1 7 8 0 0 - 1 2 0 t, 0, 1 8 0 0) \tag {13}\]

2）来袭导弹M1的运动轨迹

已知导弹M1初始位置 $P_{M1,0} = (20000,0,2000)$ ，假目标位置 $P_{\mathrm{false}} = (0,0,0)$ ，速度 $v_{M1} = 300m / s$ 。

可知导弹M1初始位置到目标位置的位移大小为:

\[\vert \vec {d} _ {M 1} \vert = \sqrt {(- 2 0 0 0 0) ^ {2} + (- 2 0 0 0) ^ {2}} = \sqrt {4 0 4 0 0 0 0 0 0} \approx 2 0 0 9 9. 7 5 m \tag {14}\]

初始位置 $P_{M1}(0)$ 代入式 (1) 可得速度向量:

\[\vec {v} _ {M 1} = (- 2 9 8. 5 1, 0, - 2 9. 8 5) \tag {15}\]

由式 (2) 可知 M1 位置函数:

\[P _ {M _ {1}} (t) = (2 0 0 0 0 - 2 9 8. 5 1 t, 0, 2 0 0 0 - 2 9. 8 5 t) \tag {16}\]

3）烟幕干扰弹及球状烟幕云团的运动轨迹

已知烟幕干扰弹投放时间 $t_{u} = 1.5s$ ，则其投放位置即无人机FY1在 $t = 1.5s$ 时的位置 $P_{\text{fall}}(1.5) = P_{\text{FY1}}(1.5) = (17800 - 120 \times 1.5, 0, 1800) = (17620, 0, 1800)$ ，已知烟幕干扰弹起爆时间 $t_{u} + \Delta t = 1.5s + 3.6s = 5.1s$ ，由式(4)知，起爆点位置即烟幕干扰弹在空中运动3.6s时的位置为：

\[\begin{array}{l} P _ {f a l l} (5. 1) = \left(1 7 6 2 0 - 1 2 0 \times 3. 6, 0, 1 8 0 0 - \frac {1}{2} \times 9. 8 \times 3. 6 ^ {2}\right) \\ = (1 7 1 8 8, 0, 1 7 3 6. 4 9 6) \tag {17} \\ \end{array}\]

又由式 (5) 知烟幕云团球心位置函数

\[P _ {C} (t) = (1 7 1 8 8, 0, 1 7 3 6. 4 9 6 - 3 \times (t - 5. 1)) \quad (5. 1 s \leq t \leq 2 5. 1 s) \tag {18}\]

Step2 迭代模拟求解

我们将从烟幕起爆时刻 $(t = 5.1s)$ 开始，到烟幕失效时刻 $(t = 5.1s + 20s = 25.1s)$ 结束，以一个微小的时间步长 $\mathrm{d}t$ 进行迭代计算。

在每一个时刻 $t$

计算该时刻导弹位置 $P_{M}(t)$ 和烟幕云团球心位置 $P_{C}(t)$ 。
导弹与烟幕球位置关系判断: 计算导弹与烟幕云团球心的距离 $d_{M_1C} = \vert P_{M_1}(t) - P_C(t)\vert $ 。如果 $d_{M_1C} > r$ , 则导弹在烟幕外部, 进入步骤3判定是否完全遮蔽; 如果 $d_{M_1C} \leq r$ , 则导弹已飞入烟幕内部, 判定为完全遮蔽, 累加 $dt$ 到总遮蔽时长, 并跳过第3步计算。
阴影锥法: 根据导弹位置 $P_{M_1}$ , 烟幕云团球心位置 $P_C(t)$ , 以及真目标最小外接球球心 $T$ 位置 $(0,200,5)$ , 由式 (7)(8)(9) 计算出半顶角 $\alpha, \beta$ , 夹角 $\theta$ , 以及 $\vert \vert \vec{M_i} T\vert \vert $ 和 $\vert \vert \vec{CT} \vert \vert $ 。若满足式 (11)(12), 则判定为完全遮蔽, 累加 $dt$ 到总遮蔽时长。
若 $t + dt \leq 25.1s$ , 进行下一轮循环判定, 否则结束循环。

循环结束后，得到有效遮蔽总时长 $= 1.3800s$ 为最终答案。

5 问题2的建模与求解

5.1 问题2模型建立

问题2是单目标连续约束优化问题，要求确定合理的无人机飞行方向、飞行速度，烟幕干扰弹投放和起爆时间，使得烟幕云团对导弹的有效遮蔽时长最大。题目中限制无人机飞行速度在 $70\sim 140m / s$ 。为了符合实际物理情况，我们将无人机飞行角度限制在 $[0,2\pi ]$ 区间内。为了使烟幕干扰弹在导弹打到假目标前有效遮蔽，我们限制烟雾干扰弹的投放时间和起爆时间在 $[0,T_{\mathrm{missile}}]$ 区间内，其中， $T_{\mathrm{missile}}$ 是导弹M1命中假目标所需时间，具体计算为：

\[T_{\text{missile}} = \frac {\sqrt {20000^{2} + 2000^{2}}}{300} = 67.0 \mathrm {s} \tag{19}\]

同时，为了保证烟幕干扰弹在触地前起爆，要求烟幕起爆点高度大于0，记无人机FY1起始位置的z坐标为 $z_{1,0}$ ，则：

\[z _ {1, 0} - \frac{1}{2} g \Delta t ^ {2} > 0 \tag{20}\]

根据求解目标和约束条件，具体建模如下：

决策变量：无人机飞行速度 $v$ ，无人机飞行方向与 $x$ 轴正方向夹角 $\phi$ ，无人机从初始位置开始飞行到投放烟幕干扰弹的时间 $t_u$ ，干扰弹从投放到起爆的滞空时间 $\Delta t$ ，即：

\[X = \left[ v, \phi , t _ {u}, \Delta t \right] \tag{21}\]

目标函数：目标是求得最大有效遮蔽时长 $T$ 。因此，目标函数是：

\[\max_ {X} f(X) = T \tag {22}\]

约束条件：决策变量必须满足以下约束条件

\[\left\{ \begin{array}{l} 7 0 \leq v _ {i} \leq 1 4 0, \text {无人机速度约束} \\ 0 \leq \theta_ {i} < 2 \pi , \text {无人机航向约束} \\ \Delta t \geq 0, t _ {u} \geq 0, \text {时间非负约束} \\ t _ {u} + \Delta t \leq T _ {\text {m i s s i l e}}, \text {烟幕有效性约束} \\ F Y _ {1 _ {z}} (0) - \frac {1}{2} g \Delta t ^ {2} > 0, \text {烟幕高度约束} \\ \Vert \vec {M _ {i}} T \Vert \geq \Vert \vec {C T} \Vert , \text {有效遮蔽物理约束} \\ \beta + \theta \leq \alpha , \text {有效遮蔽模型几何约束} \end{array} \right. \tag {23}\]

核心的“阴影锥遮蔽几何约束”本身不是一个简单的边界约束，而是内置于目标函数 $f(X)$ 的计算过程中。对于一组变量 $X$ ，我们通过问题1求解步骤2的迭代模拟方法，逐步叠加微小时间步长 $\mathrm{d}t$ ，得到遮蔽时间 $T$ ，这个过程已经包含了对物理遮蔽条件的判断。

5.2 问题2模型求解

从本问起，决策变量变多，且决策变量变化范围大，故解空间复杂，若直接根据目标函数求其解析解，则可能存在难以在有限时间里完成求解、陷入局部最优解等问题，故采用粒子群算法来实现对其求解。下面先对粒子群算法进行简要介绍。

5.2.1 粒子群算法简要介绍

粒子群算法（Particle Swarm Optimization, PSO）是一种基于群体智能的优化算法，由美国社会心理学家James Kennedy和电气工程师Russell Eberhart于1995年共同提出，现已成为计算智能领域最流行的优化算法之一。

粒子群算法模拟了群体中个体之间的社会协作和信息共享。算法将每个潜在解看作搜索空间中的一个“粒子”，整个粒子群在解空间中飞行寻找最优解。每个粒子根据自身的历史最佳位置和群体中其他粒子发现的最佳位置来调整自己的飞行方向和速度[3]。

粒子群算法的优势在于其强大的全局搜索能力和良好的收敛特性，特别适合处理高维、非线性、不可微的复杂优化问题。相比传统的基于梯度的优化方法，粒子群算法不需要目标函数的梯度信息，特别适合本题的情况。

5.2.2 利用粒子群算法的求解流程

下面我们使用粒子群优化算法解决本题，算法步骤如算法1所示。

使用粒子群算法模拟得到最大遮蔽时间是4.60秒，干扰弹需要在0秒投放，并在0.8934秒爆炸。这个解在精度上并不令人满意，使用模拟算法时可以用减小时间步长的方法增大精度，但是伴随着时间步长的减小，计算所需的资源将大大增加。所幸的是，由于我们使用的”阴影锥遮蔽模型”有足够精确的遮蔽判定条件，因此我们可以使用数值计算的方法得到更加精确的数值解。

5.2.3 数值求解步骤

定义临界函数: 我们的目标是找到真目标包围球恰好被阴影遮盖的精确时刻, 即遮蔽条件等式 $\theta (t) + \beta (t) = \alpha (t)$ 成立的时刻。我们可以定义临界函数

\[\operatorname {m a r g i n} (t) = \alpha (t) - (\theta (t) + \beta (t))\]

当 $\operatorname{margin}(t) > 0$ ，目标被遮蔽。
当 $\operatorname{margin}(t) < 0$ ，目标未被遮蔽。
当 $\operatorname{margin}(t) = 0$ ，这就是遮蔽状态变化的精确时刻。

粗略扫描定位区间: 我们仍然使用一个较粗的时间步长 $(0.1 s)$ 来扫描整个可能的遮蔽窗口。但我们的目的不是累加时间, 而是找到 $\operatorname{margin}(t)$ 符号发生变化的粗略区间 $[a, b]$ 。例如, 如果在 $t = a s$ 时 margin 为正, 而在 $t = b s$ 时变为负, 那么我们就可以确定, 遮蔽结束的精确时刻就在 $[a, b]$ 这个区间内。
数值计算求根: 使用数值求根算法来对包含临界点的区间 $[a, b]$ 内寻找 $\mathrm{margin}(t) = 0$ 的精确解。本团队采用布伦特方法求解。布伦特方法在每一步会优先尝试速度较快的割线法或逆二次差值法计算猜测根; 如果猜测根离开安全阈值, 将舍弃猜测值并退回上一步, 使用较为稳定的二分法计算数值解。

使用数值方法得到的精确的有效遮蔽时长为 $4.56897970$ 秒。此时无人机FY1的飞行角度为 $8.8122$ 度，飞行速度为70米每秒；烟幕干扰弹投放点坐标为 $(17800.00,0.00,1800.00)$ ，烟幕干扰弹起爆点坐标为 $(17861.80,9.58,1796.09)$ 。

6 问题3的建模与求解

6.1 改进阴影锥的有效遮蔽判定

问题3与问题2的主要区别是投放的烟雾干扰弹数目由1个增加为3个，除了单个球状烟幕云团的有效遮蔽，还存在两个球状烟幕云团部分重合后产生的有效遮蔽，故需改进阴影锥法。

(1) 观测平面的确立

过导弹 $M_{1}$ 与真目标中心 $T$ 的直线记为 $l_{M_{1} T}$ , 过真目标中心作与直线 $l_{M_{1} T}$ 垂直的平面即为观测平面 $\Pi$ 。

(2) 烟幕云团在观测平面上的等大投影计算

将烟幕云团 $S_{i} (i = 1,2)$ 在观测平面 $\Pi$ 上等大投影为半径为 $r = 10$ 圆形，记 $O_{i}$ 为烟幕云团 $S_{i}$ 的球心 $C_{i}$ 在 $\Pi$ 上的正投影，坐标由投影公式计算：

\[O _ {i} = C _ {i} - \frac {\vec {l} \cdot (C _ {i} - T)}{\vec {l} \cdot \vec {l}} \vec {l} (i = 1, 2) \tag {24}\]

(3) “判定小球”相关参数计算

计算两投影圆的圆心距离 $d = \Vert O_1 - O_2\Vert $ ，当 $\vert R_{O1} - R_{O2} \vert < d < R_{O1} + R_{O2}$ ，即

$0 < d < 20 \mathrm{~m}$ 时, 两圆相交, 交点连线的线段长度为:

\[L = 2 \sqrt {r ^ {2} - \left(\frac {d}{2}\right) ^ {2}} \tag {25}\]

定义“判定小球”为以 $C_1$ 和 $C_2$ 中点 $C_{\mathrm{judge}}$ 为圆心、半径 $s = \frac{L}{2} = \sqrt{r^2 - \left(\frac{d}{2}\right)^2}$ 的球体。

(4) 分层遮蔽判定逻辑

(I) 一级判定: 单烟幕云团阴影锥遮蔽判定

对每个烟幕云团 $S_{i}$ $(i = 1,2)$ , 按原来阴影锥法判定真目标是否被完全遮蔽:

(i) 计算阴影锥半顶角 $\alpha_{i} = \arcsin \left(\frac{r}{\Vert M_{1}(t)C_{i}\Vert}\right)$
(ii) 计算目标视角 $\beta = \arcsin \left(\frac{R_{\mathrm{target}}}{\Vert M_1(t)T\Vert }\right)$
(iii) 计算偏离角 $\theta_{i} = \arccos \left(\frac{(M_{1}(t)C_{i})\cdot(M_{1}(t)T)}{\Vert M_{1}(t)C_{i}\Vert \cdot\Vert M_{1}(t)T\Vert }\right)$ ;
(iv) 若满足 $\Vert M_1(t)T\Vert > \Vert M_1(t)C_i\Vert $ 且 $\theta_{i} + \beta \leq \alpha_{i}$ , 则判定 $S_{i}$ 可完全遮蔽真目标, 整体遮蔽有效。

(II) 二级判定: “判定小球” 遮蔽判定

若一级判定中 $S_{1} 、 S_{2}$ 均无法完全遮蔽真目标, 则判定真目标最小包围球 $T$ 是否被 “判定小球” $S_{\text {judge }}$ 和导弹构成的新阴影锥完全覆盖。沿用阴影锥法判定:

(i) 计算判定小球阴影锥半顶角: $\alpha_{\mathrm{judge}} = \arcsin \left(\frac{s}{\Vert M_1(t)C_{\mathrm{judge}}\Vert }\right)$ ;
(ii) 计算目标视角: $\beta = \arcsin \left(\frac{\sqrt{74}}{\Vert M_1(t)T\Vert }\right)$ ;
(iii) 计算偏离角: $\theta_{\mathrm{judge}} = \arccos \left(\frac{(M_1(t)\dot{C}{\mathrm{judge}})\cdot(M_1(t)T)}{\Vert M{1}(t)C_{\mathrm{judge}}\Vert \cdot\Vert M_1(t)T\Vert }\right)$ ;
(iv) 若满足 $ \Vert M_{1}(t)T \Vert > \Vert M_{1}(t)C_{\mathrm{judge}} \Vert $ 且 $ \theta_{judge} + \beta \leq \alpha_{judge} $ ，则判定双烟幕组合可完全遮蔽真目标，整体遮蔽有效。否则判定真目标无法被有效遮蔽。

(5) 排除三个烟雾云团重叠的情况

考虑第一个烟雾云团和第三个烟雾云团距离最近的情况，三个烟雾干扰弹间隔 1s 依次投放，并且在第三个烟雾干扰弹释放的瞬间三个烟雾干扰弹同时起爆，此时第一个云团和第三个云团的距离为 $\frac{1}{2} gt^2 = 19.6m$ ，而两个烟雾云团有重叠时球心距离需要小于 $20\mathrm{m}$ ，二者数值非常接近。为了有效遮蔽时长最大化，第一个烟雾云团和第三个烟雾云团不会在同一个水平位置起爆，因而在实际求解过程中，不存在三个烟雾云团重叠的情况。

6.2 问题3模型的建立

决策变量：本问题有8个决策变量，包括无人机飞行速度 $v$ 、飞行方向 $\theta$ 、烟幕弹的投放时刻 $t_{drop,i}$ 以及投放后的起爆延迟 $\Delta t_i$ ，其中烟幕干扰弹序号 $i\in {1,2,3}$ 。

目标函数：优化目标是使得导弹 M1 的总有效遮蔽时长最大。设 $\mathcal{T}_i$ 为第 $i$ 枚烟幕弹能够成功遮蔽导弹视线的时间点集合。由于任意一枚烟幕弹形成遮蔽即有效，总有效遮蔽时间为所有单个有效时间区间的并集。目标函数是该并集的总长度：

\[\max f (X) = \operatorname {L e n g t h} \left(\bigcup_ {i} ^ {3} \left[ t _ {s, i}, t _ {e, i} \right]\right) \tag {26}\]

约束条件：决策变量必须满足以下对所有 $i \in {1,2,3}$ 均成立的约束条件：

\[\left\{ \begin{array}{l} 7 0 \leq v \leq 1 4 0, \text {无人机速度约束} \\ 0 \leq \theta < 2 \pi , \text {无人机航向约束} \\ t _ {d r o p, i} \geq 0, \Delta t _ {i} \geq 0, \text {时间非负约束} \\ t _ {u, i + 1} - t _ {u, i} \geq 1, \text {同机投放间隔约束} \\ t _ {d r o p, i} + \Delta t _ {i} + 2 0 \leq 6 7. 0, \text {烟幕有效性约束} \\ z _ {1, 0} - 4. 9 \left(\Delta t _ {i}\right) ^ {2} - 6 0 \geq 0, \text {烟幕高度约束} \\ \text {对} \forall t \in \mathcal {T} _ {i} \text {须满足遮蔽几何条件} \end{array} \right. \tag {27}\]

6.3 问题3模型求解

根据上面的分析，我们对问题3设计粒子群优化算法，如算法2所示。

问题三的粒子群算法求解结果如图3所示，即最大遮蔽时长为6.38秒。此时无人机面向7.72度飞行，速度112.72米每秒。烟幕干扰弹1在0.00秒投出，投出后立刻起爆；烟幕干扰弹2在1.00秒投出，投出后立刻起爆；烟幕干扰弹3不发挥作用。

图3 问题三最优解示意图

7 问题4的建模与求解

7.1 问题4模型建立

决策变量：本问题有12个决策变量，包括无人机飞行速度 $v_{i}$ 、飞行方向 $\theta_{i}$ 、烟幕弹的投放时刻 $t_{drop,i}$ 以及投放后的起爆延迟 $\Delta t_{i}$ ，其中无人机序号 $i \in {1,2,3}$ 。

目标函数: 优化目标是使得导弹 M1 的总有效遮蔽时长最大。设 $\mathcal{T}_{i}$ 为第 $i$ 枚烟幕弹能够成功遮蔽导弹视线的时间点集合。由于任意一枚烟幕弹形成遮蔽即有效, 总有效遮蔽时间为所有单个有效时间区间的并集。目标函数是该并集的总长度:

\[\max f (X) = \operatorname {L e n g t h} \left(\bigcup_ {i} ^ {3} \left[ t _ {s, i}, t _ {e, i} \right]\right) \tag {28}\]

约束条件：决策变量必须满足以下对所有 $i \in {1,2,3}$ 均成立的约束条件：

\[\left\{ \begin{array}{l} 7 0 \leq v _ {i} \leq 1 4 0, \text {无人机速度约束} \\ 0 \leq \theta_ {i} < 2 \pi , \text {无人机航向约束} \\ t _ {\text {d r o p}, i} \geq 0, \Delta t _ {i} \geq 0, \text {时间非负约束} \\ t _ {\text {d r o p}, i} + \Delta t _ {i} + 2 0 \leq 6 7. 0, \text {烟幕有效性约束} \\ z _ {i, 0} - 4. 9 \left(\Delta t _ {i}\right) ^ {2} - 6 0 \geq 0, \text {烟幕高度约束} \\ \text {对} \forall t \in \mathcal {T} _ {i} \text {须满足遮蔽几何条件} \end{array} \right. \tag {29}\]

7.2 问题4模型求解

根据上面的思路，我们对问题4设计粒子群优化算法，如算法3所示。

在代码实际运行中，由于模型复杂度很大，我们需要设置一个粗略的解来初始化粒子群算法，以加快算法收敛速度、提高解的质量。在寻找初始粗略解时，可以假设无人机 $FY1$ 、无人机 $FY2$ 和无人机 $FY3$ 分别单独干扰导弹 $M1$ ，此时算法退化为问题2的情况。使用粗略解启动粒子群算法后，逐步减小时间步长、调整粒子群算法超参数、增大算法精度，最终得到模拟算法结果，此时情况如下如图4所示。

图4 问题四最优解示意图

根据图像和计算结果分析，三架无人机的有效遮蔽时间窗口相互独立、投放策略互不干涉，此时，每一架无人机对导弹 $M1$ 的干扰情况退化为问题2。因此，我们可以使用问题2的数值分析方法求得较为精确的数值解。最终计算得到，无人机FY1干扰导弹 $M1$ 的有效遮蔽时间为4.577357秒，无人机FY2干扰导弹 $M1$ 的有效遮蔽时间为3.907236秒，无人机FY3干扰导弹的有效遮蔽时间为2.847201秒。综上所述，第四问的最终结果为11.331794秒。

8 问题 5 的建模与求解

8.1 问题5模型建立

对于这个复杂的决策问题，直接使用粒子群算法寻找最优解很有可能导致计算复杂度较大、无法得到全局最优解。为了简化问题，我们考虑将问题5拆分为数个小问进行求解。计算无人机初始位置到导弹轨迹线的距离，我们不难得到如表格2所示的结果。按照就近分配的原则，将每架无人机分配给距离最近的导弹。因此，无人机FY1干扰导弹 $M1$ ，无人机FY2和无人机FY4干扰导弹 $M2$ ，无人机FY3和无人机FY5干扰导弹 $M3$ 。

在分组的基础上，我们需要对每一组规划运动和投放策略。无人机 $FY1$ 干扰导弹 $M1$ 的情况与问题3情况相同。无人机 $FY2$ 和无人机 $FY4$ 干扰导弹 $M2$ 、无人机 $FY3$ 和无人机 $FY5$ 干扰导弹 $M3$ 两种情况则可以抽象为下面的情形。

决策变量：决策变量共有16个，包括无人机FYi的飞行速度与飞行方向、第j枚烟幕弹的投放时刻与引爆延迟。

目标函数: 对于给定的一个无人机 $i$ , 我们可以计算出每一枚烟幕弹 $j$ 能够产生有效遮蔽的时间区间 $[t_{s, i, j}, t_{e, i, j}]$ 。如果某枚弹在任何时候都无法形成有效遮蔽, 则其区间为空。那么, 目标函数 $f(X)$ 就是这六个时间区间的并集的总长度。

\[\max f (X) = \operatorname {L e n g t h} \left(\bigcup_ {i, j} ^ {6} \left[ t _ {s, i, j}, t _ {e, i, j} \right]\right) \tag {30}\]

粒子群算法的目标是在决策变量构成的16维搜索空间中，使有效遮蔽时长最大化、尽量不重叠。

约束条件：

\[\left\{ \begin{array}{l} 7 0 \leq v \leq 1 4 0, \text{飞行速度约束} \\ 0 \leq \phi \leq 2 \pi , \text{飞行方向约束} \\ t_{u} \geq 0, \text{投放时间非负} \\ \Delta t \geq 0, \text{引爆延迟非负} \\ t_{u, i + 1} - t_{u, i} \geq 1, \text{同机投放间隔约束} \\ t_{u} + \Delta t \leq T_{\text {missile}}, \text{烟幕弹的起爆时刻必须在导弹命中目标之前} \\ t_{d r o p, i} + \Delta t_{i} + 2 0 \leq 6 7. 0, \text {烟幕有效性约束} \\ z_{1, 0} - 4. 9 (\Delta t_{i}) ^ {2} - 6 0 \geq 0, \text{烟幕高度约束} \\ \text{对} \forall t \in \mathcal{T}_{i} , \text {须满足遮蔽几何条件} \end{array} \right. \tag {31}\]

8.2 问题5模型求解

无人机 $FY1$ 干扰导弹 $M1$ 的过程与问题 3 完全相同，不予赘述。下面以无人机 $FY2$ 和无人机 $FY4$ 干扰导弹 $M2$ 为例设计粒子群算法，算法设计见算法 4 所示。

在代码实际运行中，由于模型复杂度很大，我们需要设置一个粗略的解来初始化粒子群算法，以加快算法收敛速度、提高解的质量。在寻找初始粗略解时，可以假设无人机 $FY2$ 和无人机 $FY4$ 分别干扰导弹 $M2$ ，此时算法退化为问题3的情况。使用粗略解启动粒子群算法后，逐步减小时间步长、调整粒子群算法超参数、增大算法精度，最终得到计算结果。

无人机FY1干扰导弹M1的最大有效遮蔽时长为6.38秒，无人机FY2和无人机FY4干扰导弹M2的最大有效遮蔽时长为11.58秒，无人机FY3和无人机FY5干扰导弹M3的最大有效遮蔽时长为6.6秒。即问题5的最终结果为24.56秒，具体投放策略详见result3.xlsx附件。

9 模型误差分析

9.1 阴影锥模型的误差

在对真目标圆柱体建模时，有很多种不同的方法来表现烟雾云团“遮蔽”真目标。在这一节，我们以问题一为背景，通过实际实现不同的建模方法，证明阴影锥模型在对“遮蔽”建模上的高精度与低计算量。

建模方法一：将真目标圆柱简化为一个点

由于导弹到圆柱体的距离为几千米，相对而言圆柱体的直径可以忽略不计。因此将整个圆柱体简化为中心点，如果中心点被遮挡，则认为整个圆柱被遮挡。这个方法计算量极小，但是当导弹飞行时间增加后会产生严重的误差。

建模方法二：在真目标圆柱上取16个采样点

用目标 $T$ 表面上一组有限的采样点集合 ${\vec{p}_1,\vec{p}_2,\dots ,\vec{p}_n}$ 来近似其几何轮廓。只有当所有采样点都被遮挡时，才判定整个目标被遮挡。对于每一个采样点 $\vec{p_i}$ ，其对应的视线线段 $\overline{\vec{P_M p_i}}$ 都必须与烟幕球体 $S$ 有交集

数学上，遮蔽条件为：

\[\forall \vec {p}_{i} \in \{\vec {p}_{1}, \dots , \vec {p}_{n} \} \subset T, \quad \exists t_{i} \in [ 0, 1 ] \quad \mathrm {s . t .} \quad \Vert (\vec {P}_{M} + t_{i} (\vec {p}_{i} - \vec{P}_{M})) - \vec{P}_{S} \Vert \leq R_{S}\]

可以预见，采样法精度将大于中心点法；然而采样法的精度和计算成本与采样点数 $n$ 直接相关，在本题的背景下，如果使用采样法建模 + 粒子群算法，计算量将极大上升。

建模方法三：阴影锥 + 目标包围球法

本方法的具体原理已在问题2的建模中说明。阴影锥法通过判断锥角的等式关系与物体距离，能够快速计算出精度较高的结果。但是由于使用了目标包围球近似，存在微小的系统误差，计算得到的结果将偏小。

建模方法四：高密度点采样法

此方法是建模方法二的极限情况，在本次实验中我在圆柱上取了720个点来获得一个接近连续模型的、较精确的数值解。这个方法作为其他方法的真实值衡量误差

根据表格3，我们的阴影锥 + 目标包围球方法在不损失太多精确度的情况下，大大提升了计算速度、降低了计算复杂度。在本题粒子群算法的计算背景下，阴影锥方法不仅在数学上更优雅，也在计算效率和物理真实性之间取得了最佳平衡。

9.2 粒子群算法误差分析

粒子群算法的主要误差来源是“寻优误差”，即算法找到的解可能并非理论上的全局最优解。在代码运行的不同时刻，我会多次调整超参数以确保算法不陷入局部最优解。惯性权重过大、认知学习因子过大等都会造成探索性强但局部优化能力弱，影响算法的求解能力。

此外，粒子群算法是一种随机搜索算法，其初始化和速度更新都包含随机数。这意味着，即使参数和问题完全相同，两次独立运行也可能得到不完全一样的最优解。针对这种现象，我们使用了有初始点的粒子群算法，虽然能够加快代码速度、提高解的质量，但是也有可能使整个粒子群过快地趋向了某个局部最优解，导致种群多样性丧失。

在算法选择上，遗传算法的全局搜索能力强，不易早熟，但后期收敛速度可能较慢；而差分进化算法的鲁棒性强，收敛速度和全局搜索能力之间有很好的平衡，但是对参数的敏感性相对较低；粒子群算法通常收敛速度较快，但容易陷入局部最优；改进后的粒子群算法较好的综合了粒子群算法计算速度快的特点，同时探索能力强，能够较快地搜索到全局最优解。

10 模型评价

10.1 模型的优点

优点1：计算结果精度较高。使用粒子群算法得到精度较高的模拟解，使用数值计算方法得到精度极高的有效遮挡时间长度，答案精度很高。
优点2：计算速度快，能使用极小的计算资源得到较好的结果。
优点3：针对纯粒子群方法实验随机性强的缺陷，先讨论简化问题下的特解，再将特殊解作为初始粒子进入粒子群算法，以加快求解速度、提高求解质量。例如，问题5中对单架无人机分别使用粒子群算法，将特殊结果作为整体粒子群算法的初值，答案精度较高，极大缩短了运算时间。

10.2 模型的缺点

缺点 1: 代码随机性较强, 实验复现困难。由于粒子群算法是一种随机搜索算法, 即使参数和问题完全相同, 两次独立运行也可能得到不完全一样的最优解。这对建模问产生不利影响。
缺点2：建模复杂度不足。问题五的建模中，无人机与导弹的分组仅仅依靠初始距离，同时也忽略了单个烟雾云团对多个导弹的干扰、多个烟雾云团对多个导弹的干扰等较复杂的情况。

11 模型推广与改进

阴影锥 + 目标包围球的方法可以运用在计算机图形学与游戏开发中，为了提高渲染效率，游戏引擎可以剔除掉被其他物体遮挡而看不到的物体，此时可以使用该方法快速判断是否需要渲染光影背后的复杂模型；在机器人自动驾驶中，可以快速判断车辆背后的视觉盲区有多大、是否有可能有行人隐藏在视觉盲区中。当然，阴影锥 + 目标包围球方法也存在较大的误差，当目标的形状非常不规则时，包围球会包含大量无效空间，导致精度下降；在现实中，观察者或光源存在有体积的，这会在本影周围产生一个模糊的、部分遮蔽的区域，即“半影”。可以通过这些方法改进阴影锥 + 目标包围球方法，应用在更多精度要求不是特别高、要求快速反应判断的场景中。

参考文献

[1] 陈浩, 高欣宝, 李天鹏, 等. 国外烟幕干扰弹发展及关键技术研究[J]. 飞航导弹, 2017 (12):71-74.
[2] CROW F C. Shadow algorithms for computer graphics[J]. ACM SIGGRAPH Computer Graphics, 1977, 11(2):242-248.
[3] 刘刚, 彭春华, 相龙阳, 等. 采用改进型多目标粒子群算法的电力系统环境经济调度 [J]. 电网技术, 2011, 35(07):139-144.
[4] QWEN. Qwen3-235b-a22b-2507[R]. 杭州: 阿里巴巴集团控股有限公司, 2025.9.6.
[5] GEMINI. Gemini 2.5 pro[R]. USA: Google DeepMind, 2025.9.6.