低空突防是无人机（unmanned aerial vehicle, UAV）作战中的重要环节^[1]，其目标是在满足各项约束的基础上，设计算法规划出从起点到终点能够有效避开威胁源的低能耗路线. 单无人机突防能力略显不足，因此多无人机协同规划逐渐成为无人机集群化、自主化和智能化的发展关键和研究热点^[2]. 由于作战环境形势交错、复杂性强，且多无人机在航迹规划时受到众多条件的约束，亟须建立更精确的模型并设计更高效的算法. 元启发式算法被认为是处理单目标和多目标问题的标准优化方法^[3]，这些方法在复杂和动态环境（包括战场、救援任务和多无人机路径优化等场景）中表现出色. 元启发式算法大致分为4类：进化算法^[4]、基于物理的算法^[5]、基于群体的算法^[6]和混合元启发算法^[7]. 其中基于群体的算法也称群体智能（swarm intelligence, SI）^[8]算法，它不仅克服了进化算法收敛慢和基于物理的算法受模型限制的缺陷，而且计算成本适中. 此外，该方法强调合作和信息共享处理集体智能问题，对初始条件的敏感性较低，适用于三维多无人机低空突防. 群体智能算法已被应用于多无人机航迹规划. Eun等^[9]在二维地图上利用遗传算法讨论多无人机组的任务分配和轨迹规划问题，但未充分考虑三维地形和威胁源约束的影响. 杜云等^[10]考虑无人机性能、威胁源、时空协同等约束，利用改进的粒子群算法求解最优路径集合，但利用单目标代价函数寻优时，权重设置受主观因素影响较大. 闫少强等^[11]提出改进麻雀搜索算法（game predatory sparrow search algorithm, GPSSA），在复杂三维地形中求解单无人机及多无人机航迹规划中都取得了较好的规划效果，但协同目标函数仍依赖单目标群智能算法的寻优性能，影响了算法效率. 袁梦顺等^[12]将各项约束和代价转换为非支配排序遗传算法-Ⅲ(non-dominated sorting genetic algorithm-Ⅲ, NSGA-Ⅲ) 的数值指标，避免了单目标问题加权的主观性，但算法应用的建模场景复杂程度不高. 上述研究虽然能够实现多无人机协同航迹规划，但优化效果和迭代收敛速度仍有提升空间.

秃鹰搜索算法(bald eagle search，BES)^[13]是新型元启发式算法，具有收敛速度快，适应性强，模型易修改等特点. 与同类型的如麻雀算法(sparrow search algorithm，SSA)、鹈鹕算法(pelican optimization algorithm，POA)、水母算法(jellyfish search，JS)和蜣螂算法(dung beetle optimizer，DBO)相比，BES的性能更优，且在极坐标系中进行搜索的寻优方式独特，更适用于高效搜索无人机构型空间. 多无人机航迹规划问题通常涉及多个冲突的目标，使用单目标群智能算法只能解决单个目标或通过加权方式平衡多个目标，求解精度欠缺. 多目标算法能够在优化过程中同时优化多个目标，找到一系列平衡解，提供多样化决策选择. BES本身搜索能力不足，寻优过程易趋于中心点，选取局部最优航迹时缺陷较大. 为了使多无人机在航迹规划中实现快速寻优，获取到具有高精度的较优路径，本研究提出改进的多目标秃鹰搜索算法(improved multi-objective bald eagle search algorithm，IMBES)，应用于多无人机协同低空突防.

假设静态地形环境已知，根据已知参数构建三维山地模型、威胁源模型、单无人机物理约束模型、多无人机协同约束模型以及路径平滑度约束模型，结合各约束确立多目标代价函数.

建模时地形约束有2种，分别是地势低平的平原和地势较高的山丘^[14-15]，其中山丘对地形代价影响较大. 构建山丘三维模型为

(1)$ H(x,y) = h \cdot \exp {\left( { - \frac{{{{(x - {x_0})}^2}}}{{{\omega _1}}} - \frac{{{{(y - {y_0})}^2}}}{{{\omega _2}}}} \right)^2}{\text{.}} $

式中：(x, y)为三维空间点的坐标，(x₀, y₀)为山丘的中心坐标，h为山丘高度，ω₁、ω₂均为山丘陡峭度. 无人机飞行过程中，离地面高度高于0.1 km即视作安全. 航迹点地形威胁代价函数为

(2)$ f_{_{H,j}}^i = \left\{ \begin{gathered} K_{_H}^i{\text{,}}\quad h_j^i - H_j^i \leqslant 0{\text{;}} \\ K{_{_H}^{i'}}{\text{,}}\quad 0 < h_j^i - H_j^i \leqslant 0.1{\text{;}} \\ 0,\;\;\;\quad h_j^i - H_j^i > 0.1. \\ \end{gathered} \right. $

式中：$f_{_{H,j}}^i $、$h_j^i $、$H_j^i $、$K_H^i $及$K{_{_H}^{i'}} $分别为第i架UAV第j个航迹点对应的地形代价、海拔高度、地形高度和2个相异的地形威胁系数. 为了更大程度地接近实际，低空突防须考虑高程代价. 各航迹点的高程代价为

(3)$ f_{_{G,j}}^i = h_{_j}^i. $

限定无人机飞行边界以保证可控性，根据地形情况设定飞行水平范围及最大飞行高度.

在建模无人机突防过程中，雷达探测、防空火炮、地空导弹和禁飞区等被认为是主要的威胁源. 为了简化建模，根据威胁的性质算出各威胁源的作用半径和作用空间形状，以等效地形代替^[16]. 假设距离威胁源中心越近，威胁代价越大；反之，则越低. 将威胁源等效为数学模型，定义每个航迹点在威胁区域的计算代价为

(4)$ f{_{P,j}^{m,i}} = \left\{ \begin{array}{*{20}{l}} 0,& | {(x,y)}_j^i - {C_n^m} | \geqslant {R^m_{\max,n}}; \\ K_{_{{\mathrm{thr}}}}^m| {R_{\max n}^m - | {(x,y)} _j^i - {C_n^m} |} |,& 其他. \\ \end{array} \right. $

式中：$K_{{\mathrm{thr}}}^m $为第m类威胁源的修正系数，$C_n^m $、$R_{{\mathrm{max}},n}^m $分别为第m类的第n个威胁源的中心及其最大作用半径，$(x,y)_j^i $、$ f{_{P,j}^{m,i}} $为第i架UAV第j个航迹点对应的笛卡尔坐标和威胁源代价.

单无人机在飞行过程中受自身的物理约束，主要包括最大拐弯角α、最大下滑及爬升角β、续航时间t_d以及发动机功率W约束^[17]，表达式分别为

(5)$ J_{{\alpha},j}^i = \left\{ {\begin{array}{*{20}{l}} {0,}&{\alpha \leqslant \alpha _{_{\max }}^i;} \\ {{K_\alpha },}&{\alpha > \alpha _{_{\max }}^i.} \end{array}} \right. $

(6)$ {J}_{\beta,j}^{i}=\left\{\begin{array}{*{20}{l}}0\text{，}&\beta \leqslant {\beta }_{{}_{\mathrm{max}}}^{i};\\ {K}_{\beta },& \beta > {\beta }_{{}_{\mathrm{max}}}^{i}.\end{array}\right. $

(7)$ J_{t_{{{\mathrm{d}}}},j}^i = \left\{ {\begin{array}{*{20}{l}} {0,}&{{t_{\mathrm{d}}} \leqslant {t}_{_{\max }}^i;} \\ {{K_{{t_{\rm{d}}}}},}&{{t_{\mathrm{d}}} > {t}_{_{_{\max }}}^i.} \end{array}} \right. $

(8)$ J_{{W},j}^i = \left\{ {\begin{array}{*{20}{l}} {0,}&{W \leqslant W_{_{_{\max }}}^i;} \\ {{K_W},}&{W > W_{_{_{\max }}}^i.} \end{array}} \right. $

式中：K_α、K_β、$K_{{t_{\rm{d}}}} $、K_W分别为拐弯角、俯仰角、续航时间和无人机发动机功率的威胁系数；$ \alpha _{_{\max }}^i $、$ \beta _{_{\max }}^i $、$ {t}_{_{\max }}^i $和$W_{_{\max }}^i$分别为第i架UAV对应的最大拐弯角、俯仰角、续航时间和无人机发动机功率；$ J_{{\alpha},j}^i $、$ J_{\beta ,j}^i $、$ J_{t_{{{\mathrm{d}}}},j}^i $和$J_{{W},j}^i$分别为第i架UAV第j个航迹点对应的不同约束的代价函数. 第i架UAV第j个航迹点自身约束代价函数为

(9)$ f_{J,j}^i = J_{{\alpha},j}^i+J_{{\beta},j}^i+J_{t_{{{\mathrm{d}}}},j}^i+J_{{W},j}^i. $

路径平滑是生成路径可行的必要条件，通过计算转弯速率和爬坡速率得到^[18]. 如图1所示，转弯角$\phi_j^i $为投影在水平面oxy的2个连续路径段$ \rightharpoonup {P{{_j^{i'}}}P{{_{j+1}^{i'}}}} $和$ \rightharpoonup {P{{_{j+1}^{i'}}}P{{_{j+2}^{i'}}}} $之间的夹角. 设k为z轴方向的单位向量，则投影向量计算式为

(10)$ \rightharpoonup{P{{_j^i}^{'}}P{{_{j+1}^{i'}}}} = {\boldsymbol{k}} \times (\rightharpoonup {P_j^iP_{j+1}^i} \times {\boldsymbol{k}}). $

转弯角为

(11)$ {\phi }_{j}^{i}=\mathrm{arctan}\left(\frac{\left\| \rightharpoonup{{P}_{j}^{i'}{P}_{j+1}^{i'}}\times \rightharpoonup{{P}_{j+1}^{i'}{P}_{j+2}^{i'}}\right\| }{\rightharpoonup{{P}_{j}^{i'}{P}_{j+1}^{i'}}·\rightharpoonup{{P}_{j+1}^{i'}{P}_{j+2}^{i'}}}\right). $

爬升角$\varPsi_j^i $是路径段$ \rightharpoonup {P_j^iP_{j+1}^i} $与该路径段在水平面上的投影$ \rightharpoonup {P{{_j^i}^{'}}P{{_{j+1}^{i'}}}} $的夹角，计算式为

(12)$ \psi _j^i = \arctan \left( {({{z_{j+1}^i - z_j^i}})\Bigg/{{\left\| {\rightharpoonup {P{{_j^i}^{'}}P{{_{j+1}^{i'}}}} } \right\|}}} \right). $

综上，第i架UAV路径平滑度代价为

(13)$ f_{_s}^i = {\lambda _1}\sum\limits_{j = 1}^{n - 2} {\phi _{_j}^i+{\lambda _2}} \sum\limits_{j = 1}^{n - 1} {\left| {\psi _j^i - \psi _{j - 1}^i} \right|} . $

式中：λ₁、λ₂分别为转弯角和爬升角的惩罚系数.

多无人机系统在执行协同打击作战任务时，须在保证同时到达同一目标点的前提下，尽可能实现各无人机的飞行路径足够分散，即实现“分散突防、集中打击”的作战效果，从战术角度保证以最小代价实现最大打击效果^[19]. 设时间和空间协同约束系数$C_{i,j}^{\mathrm{t}} $和$C_{i,j}^{\mathrm{d}} $分别为

(14)$ C_{i,j}^{\text{t}} = \left\{ \begin{gathered} P/\left| {{t_i} - {t_j}} \right|,\max \left[ {{t_{i,\min }},{t_{j,\min }}} \right] < \min \left[ {{t_{i,\max }},{t_{j,\max }}} \right]; \\ 0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;其他. \\ \end{gathered} \right. $

(15)$ C_{i,j}^{\text{d}} = \left\{ \begin{gathered} P'/\sum\limits_{k = 1}^n {{d_{k}},} \quad {d_{k}} < 0.3{\text{ km}}; \\ 0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\quad 其他. \\ \end{gathered} \right. $

式中：P、P^'均为足够大的常数，本研究均取为100；|t_i－t_j|为无人机i和无人机j到达目标点的时间差；$C_{i,j}^{\mathrm{t}} $为无人机之间的时间协同系数；时间协同条件为max[t_i,min, t_j,min]<min[t_i,max, t_j,max]，即无人机抵达时间范围存在交叉；d_k为2架无人机第k个航迹点的欧氏距离，n为航迹点数. 假设有M架无人机，则第i架无人机与第j架无人机的协同时空约束代价函数为

(16)$ {f_{\mathrm{c}}} = 1/\sum\limits_{j = 1}^M {\left( {C_{i,j}^{\text{t}}+C_{i,j}^{\text{d}}} \right)} . $

多无人机协同航迹规划问题，既包含单个无人机的路径规划，又涉及各机之间的时空协同性. 选取各无人机的总路径长度l_i和综合威胁代价f_i这2个相互冲突的目标构建适应于多无人机协同的路径规划的多目标优化模型：

(17)$ \min {\left[ {{l_i},{f_i}} \right]^{\text{T}}};\;i = 1,2, \cdots ,M. $

(18)$ {l_i} = \sum\limits_{k = 1}^n {{d_{{k-1,k}}}} . $

(19)$ {f_i} = \sum\limits_{j = 1}^n {({\sigma _1}f_{_{J,j}}^i} +{\sigma _2}f_{_{H,j}}^i+{\sigma _3}f_{_{G,j}}^i+{\sigma _4}f_{_{P,j}}^{m,i})+{\sigma _5}{f_s}+{\sigma _6}{f_{{\mathrm{c}}}}. $

式中：d_k−1,k为路径点k-1和k的欧氏距离； 在实际应用中，须设置权重以满足需求σ₁、σ₂、σ₃、σ₄、σ₅、σ₆均为权重系数.

结合式(17)将航迹规划问题转化为多目标优化问题.

设定迭代中每条更新路径由种群中的秃鹰X_i={x_i1, x_i2, ···, x_id}构成. 在优化过程中，利用球面矢量的幅度、仰角和方位角分量与无人机的速度、转弯角和爬升角之间的相互关系来获取更安全、高质量的解决方案^[18]. 用球坐标系代替笛卡尔坐标表示秃鹰位置，因为所有路径的起点和终点固定，所以不包含在优化序列中，路径由n个航迹点组成，选用维度为3n的向量序列进行表示. 这些向量在球坐标系中表示为3个分量，分别为航迹距离增量ρ、仰角Ψ∈(−π/2,π/2)和方位角Φ∈(−π,π). 飞行路径Ω_i用超球面向量表示为

(20)$ \begin{split} {\boldsymbol\varOmega _i} =& \left[{{\rho _{i1}},{\psi _{i1}},{\phi _{i1}},{\rho _{i2}},{\psi _{i2}},{\phi _{i2}}, \cdots ,{\rho _{iN}},{\psi _{iN}},{\phi _{iN}}} \right],\\& N = n - 2.\end{split} $

航迹点P_i,j=(x_i,j，y_i,j，z_i,j)须转化为笛卡尔系，转化式为

(21)$ \left. \begin{gathered} {x_{i,j}} = {x_{i,j - 1}}+{\rho _{ij}}\sin {\psi _{ij}}\cos {\phi _{ij}}{\text{,}} \\ {y_{i,j}} = {y_{i,j - 1}}+{\rho _{ij}}\sin {\psi _{ij}}\cos {\phi _{ij}}{\text{,}} \\ {z_{i,j}} = {z_{i,j - 1}}+{\rho _{ij}}\cos {\psi _{ij}}. \\ \end{gathered} \right\} $

综上，在球坐标系中秃鹰根据式(17)中的代价函数，按照更新公式迭代，逐渐收敛到最优Pareto前沿，同时将秃鹰种群从角空间映射到笛卡尔坐标空间，得到笛卡尔坐标系下的完整路径.

群体智能须维护一组解并利用启发式规则搜索最优解. 初始种群是进化过程中的起始猜想，影响SI算法在寻找全局最优解时的起点，影响种群的收敛速度和最终解的精度. SI算法多采用伪随机数生成器(pseudorandom number generators，PRNGs)生成初始种群，其原理是通过PRNGs生成一组服从均匀分布的随机数，进而覆盖搜索空间中可能的区域. PRNGs的优点是算法简单，适用范围广. 混沌数生成器(chaotic number generators, CNGs)是基于混沌技术的随机数生成器，利用一维混沌映射，指定随机初值不断迭代，生成一系列连续的点，产生混乱的种群. 相比PRNGs，CNGs能够在种群多样性、成功率以及收敛性等方面提高群体智能算法性能. 混沌序列行为对初值和参数过于敏感，针对不同问题须选择不同的初值和参数. 同时，由于吸引子的存在，混沌序列可能被吸引到几个固定的点，影响算法效果.

本研究设计耦合混沌映射，如图2所示为各混沌映射分岔图，横轴为参数值，纵轴为映射的状态或输出值. 由Chebyshev映射^[20]和Infinite Fold^[21]混沌映射的分岔图可以看出，当2个混沌映射中的控制参数分别在大于1.0和1.8时，逐渐出现分岔现象，因此采用耦合参数c级联这2个种子映射，对非线性映射的混沌片段进行拼接，通过交叉迭代的方式从一维扩展到多维，其中二维情况的表达式为

(22)$ \left. \begin{gathered} {x_{i+1}} = \sin \;({\text{π }}(\cos \;(3 - c){\cos ^{ - 1}}({x_i})+\sin \;(c/{y_i})), \\ {y_{i+1}} = \sin \;({\text{π }}(\cos \;(3 - c){\cos ^{ - 1}}({y_i})+\sin \;(c/{x_i})). \\ \end{gathered} \right\} $

式中：c=0.8，此时3−c＞1.8，2个种子映射均出现混沌现象；i为迭代次数. 如图3所示，改进后的耦合映射分岔能够减少映射空间的局部偏置，生成更均匀的种群空间，提高算法的搜索性能. 将耦合混沌映射初始化种群的BES (CO-MOBES)与原始多目标BES (MOBES) 在ZDT1测试函数上进行对比实验，结果如图4所示. 图中，f₁为测试函数的目标函数1，f₂为测试函数的目标函数2. 由图可知，改进算法有效提高了初始化种群质量，pareto前沿覆盖率较高.

多目标秃鹰搜索算法在选择和搜索阶段完成后，引导整个秃鹰种群向固定方向寻找非支配解. 由于该机制依赖前期充分探索解空间的程度，导致秃鹰个体盲目跟随最优秃鹰并逼近最优位置附近的非支配解，不进行更优解探索；若此时赋予部分秃鹰个体自主探索的能力，将会提高整体开发能力. 在此阶段引入“侦察鹰”，侦察鹰搜索机制通过自适应高斯游走策略更新个体位置，逐步调整搜索方向和步长ϛ，以自适应的方式在搜索空间中进行探索求解. 设计组内均值引导，动态平衡算法探索与开发能力. 捕获阶段更新公式为

(23)$ \varsigma = = {\text{exp}}\left( {\frac{{2.5 {\mathrm{rand}} - 1}}{{{{\left( {1 - {I_{{\mathrm{max}}}}} \right)}^2}}} \ln i} \right). $

(24)$ {p}_{i,\text{n}}=\left\{ \begin{array}{*{20}{l}}G\text{(}{p}_{i,\text{m}},\varsigma \text{)}+\text{(}{r}_{1}{p}_{i,\text{m}}-{r}_{2}{p}_{i}\text{),}& {p}_{i}=\dfrac{1}{3}p;\\{\mathrm{rand}} {p}_{\text{b}}+{\delta }_{x}+{\delta }_{y}\text{,}& {p}_{i}=\dfrac{2}{3}p.\end{array} \right.$

式中：r₁、r₂均为服从[0,1.0]均匀分布的随机数，rand为随机函数，随机生成[0,1.0]的随机数， I_max为最大迭代次数，δ_x、δ_y分别为极坐标中秃鹰的变化值，p_b为当前秃鹰最好的搜索位置；p_i,m为第i次迭代秃鹰更新位置的平均值，p为种群数量，p_i,n为第i次迭代秃鹰更新后的位置. 在迭代初期, 高斯分布具有较大的步长ϛ以保证算法具有较强的探索能力；后期ϛ逐渐减小，提高了算法的局部开发能力. 为了削弱固定方向引导效应，通过组内均值个体引导，在其附近产生新的个体，增强侦查能力. 将基于“侦察鹰”的自适应高斯游走策略的BES (Inves-BES) 与MOBES在ZDT1测试函数上进行对比实验，改进前后代价变化曲线对比如图5所示. 由图可知，改进算法寻优性能进一步提升.

相比现有的多目标排序方法（如基于网格法的非支配排序），快速非支配排序策略只与已经分配到前沿的解进行比较，无须进行所有解比较，减少了计算量. 该方法的排序步骤如下：1)得到目标函数适应度值的矩阵，依次遍历判断支配关系，取其中的非支配解集，构成rank1，即得到首个Pareto前沿；2)将其他解与rank1比较，删去其他解中由rank1所有个体支配的次数，若此时该解为非支配解，那么该解位于rank2，据此逐一判断得到rank2；3)重复步骤1)、2)，将后面某个解与已经得到分配前沿的解进行比较，即可确定rank3、rank4、···、rankn.

通过改进策略，算法不仅可以快速找到最优解，而且可以在整个解空间中保持多样性，进而在应用于多无人机协同航迹规划时具有优越的性能和效果. 将引入快速非支配排序策略的BES(NSBES)与MOBES在ZDT1测试函数上进行对比实验，结果如图6所示. 可以看出，非支配排序得到的Pareto前沿解精度更高.

在无人机实际飞行中，通过三次样条平滑路径寻优航迹点. 设最终无人机避障路径L={l₁, l₂, ···, l_d}，在起点与终点之间进行三次样条插值，最终形成一条平滑航迹曲线. 针对算法运行中出现路径点过于集中的问题，前期在最优个体样条插值的平滑路径点上随机选点，使得路径点尽可能随机分布，算法流程如下. 1）种群初始化：构建三维环境模型、威胁源模型、无人机物理约束模型以及路径平滑度约束模型，确定航迹点数目与种群数目，根据式(22)利用耦合混沌映射初始化秃鹰种群，计算多目标代价函数，进行非支配排序. 2）迭代更新：根据秃鹰位置更新公式进行搜索. 3）边界限定：根据已限定的三维空间范围检测是否边界溢出，若溢出，则用该航迹点溢出维度的边界值代替. 4）三次样条插值随机取点获取平滑路径. 5）判断是否达到最大迭代次数，若未达到，则跳回步骤2；否则结束.

实验硬件仿真平台使用处理器为12th Gen Intel(R) Core(TM) i7-12700H 2.30 GHz，内存为16.0 GB的计算机，软件平台为Matlab R2020b. 设置种群规模为500，迭代次数为200. 在标准测试函数集上将IMBES与基于网格法的多目标麻雀算法(MOSSA)^[22]、多目标鹈鹕算法(MOPOA)^[23]、多目标灰狼算法(MOGWO)^[24]、基于非支配排序的多目标蜣螂算法(NSDBO)^[25]以及NSGA-Ⅱ^[26]等主流算法的求解性能进行对比测试.各测试函数具体描述如表1所示，其中N_OF为目标函数个数.

测试过程中，各测试函数分别运行30次，得到的pareto前沿箱线图如图7所示，其中R为测试结果分布范围. 可以看出，IMBES在解决各种多目标函数方面表现出色. 对于测试函数ZDT1，算法MOSSA、MOGWO和NSDBO都表现出相对较高的中位数，MOPOA在中位性能上处于中等水平. 算法NSGA-Ⅱ、MOBES和IMBES的中位数相对较低，且离散度低，表明性能相对稳定. IMBES具有最低的中位数和最佳性能. 对于Kursawe测试函数，NSGA-Ⅱ表现出总体稳定的性能，离群值较少，中位数属于中等范围. IMBES具有最低的中位数，表现出更强的鲁棒性. 对于三目标函数，所提算法收敛性和覆盖率仍占优；根据没有免费午餐定理(no free lunch, NFL) ^[27]，不能期望算法在所有测试函数都有最好的表现，因此对于DTLZ6这样较为复杂的不连续多模态测试函数，IMBES未呈现明显优势. 但IMBES整体寻优效果较优，说明IMBES在复杂交错环境中的寻优能力较强；MOBES相较于IMBES在各个测试函数均呈现出显著优势. 6种函数的测试结果证明了所提算法的优越性.

设M为目标的个数，N为种群大小，D为问题维度. 在原始BES中，3个阶段位置更新的总时间复杂度为O(ND)×3，本研究设计的耦合混沌映射初始化的时间复杂度为O((N lg N)+1)，基于“侦察鹰”的自适应高斯游走策略时间复杂度为O(1/3N)，引入的快速非支配排序时间复杂度为O($MN\sqrt {{N}} $)，拥挤度距离计算时间复杂度为O(N²). IMBES的时间复杂度为O(ND)×3+O((N lg N)+1)+O(1/3N)+ O($MN \sqrt {{N}} $)+ O(N²). 总体而言，算法的时间复杂度主要受到非支配排序速度影响.

如表2所示为不同航迹点数目n_al对各群体智能优化算法的影响，种群数目和迭代次数分别设为500和200，其中l_t为总路径长度，f_t为总威胁代价，t为耗时. 由表可知，1)增加航迹点数目提升了路径规划精度，增加了决策变量，算法因此能够更好地探索搜索空间，找到更优方案. 2)增加航迹点数目导致计算复杂度的增加，耗时变长，因此在确定航迹点数目时应权衡计算复杂度和路径规划的精度. 3)须综合考虑路径规划的精度、计算复杂度和无人机任务需求等因素来确定航迹点数目. 实验中考虑上述因素，航迹点数目设置为16.

为了验证所提模型合理性和IMBES的可行性，在仿真实验中将IMBES与各算法进行比较. 设置种群规模为500，迭代次数为200，航迹点数目为16. 各算法运行30次，其他主要参数设置如表3所示. 硬件和软件平台与第3节所述保持一致.

建立双模型增加算法的可信度，实验如下.

模型1：选取2架UAV进行仿真实验，建立100 km×100 km×5 km的模拟地形区域，实验参数如表4所示. 第一架无人机记作UAV1，航线起点为(15,20,2)，终点为(80,80,2). UAV2的航线起点为(40,10,1)，终点为(90,55,1). 2架无人机型号相同. 威胁源数量为5个，其中雷达威胁为2个，作用中心分别为(30,40)和(70,60)，作用半径为10 km，作用效果等效为螺旋线. 禁飞区、地空导弹和防空火炮数量均为1个，作用中心分别为(35,20)、(65,75)和(50,60)，作用半径分别为7、8和5 km，作用效果等效为圆柱体. 如图8所示为各算法航迹规划结果，各算法求解质量比较如表5所示，l_m为平均路径长度，f_m为平均威胁代价，l_s为最短路径长度，f_s为最小威胁代价. NSGA-Ⅱ规划出的路径较短，但多次穿过威胁源，威胁代价极高，规划效果差；其他算法均能规划出满足约束的路径. 由表5可知，所提算法稳定性更高，无论是航程代价还是威胁代价，均得到最优的规划效果. 由图8可以看出，所提算法规划的路径中，UAV1和UAV2飞行路径没有发生交叉和重合，这说明所提算法在对无人机路径协同规划问题进行寻优时实现了无人机之间的空间协同，算法的寻优效果能够让飞行路径更短，威胁代价更小，耗时最短.

模型2：测试场景基于从激光雷达传感器获得的真实数字高程模型（DEM）地图^[28]. 在此场景下，设置5个威胁. 模型2的多架无人机协同规划轨迹结果如图9所示. 如表6所示为每种算法的解质量的比较，度量值表示实验中所有无人机的求和结果，其中c_s为最小威胁代价. 实验增加基于改进多目标飞蛾-火焰优化算法(LTG-NSMFO)^[29] 、多目标水母搜索算法(MOJS)^[30]、NSGA-III^[31]作为对比算法. 由图可知，所提算法规划了最优路径并获得最优帕累托前沿. 由表6可知，所提算法在轨迹规划中表现出更高的稳定性以及最优的耗时性能，可以达到最优的规划效果. 与MOSSA、MOPOA、MOJS、NSDBO、NSGA-Ⅲ和LTG-NSMFO相比，所提算法将平均路径长度缩短了1.10~35.80 km，将平均威胁成本降低了0.72~302.49；结果表明，所提算法可以在复杂地形环境中提供可靠、高效的目标规划解决方案.

如图10所示为2架无人机在每种算法规划路径的不同轨迹点的高度比较曲线. 可以看出，所提算法的高度曲线最低，MOPOA高度变化最大且整体高程代价高；所提算法的爬升角和滑动角成本最低. MOPOA和MOJS在不同飞行路径点之间的高度有显著变化，这对无人机的性能有要求. 与MOSSA、MOPOA、MOJS、NSDBO、NSGA-Ⅲ和LTG-NSMFO相比，所提算法的无人机1平均飞行高度下降了2.57~22.68 km，无人机2的平均飞行海拔下降了1.24~10.41 km.

本研究提出改进多目标秃鹰搜索算法并应用于多无人机协同低空突防，解决了三维空间中多无人机协同航迹规划及自主避障问题. 定性和定量实验结果表明，所提算法实现了快速优化，能够为多无人机自主避障协同路径规划提供精确的路径. 此外，所设计的模型能够解决地形和威胁源的影响，便于深入探索复杂环境中的算法性能. 在未来的研究中，将测试算法在城市地区的性能，评估算法在不同测试环境和场景中的表现.