技术领域
[0001] 本
发明是一种求解无人机约束航路规划的混合烟花粒子群协同
算法，属于智能
机器人导航与控制以及智能优化技术领域。
背景技术
[0002] 路径规划是机器人导航技术的一个重要研究方向。而无人机航路规划是在考虑了航程和飞行环境等约束的综合情况下，在适当的时间内计算出无人机从起始点到目标点最优或次最优的飞行轨迹，使无人机能回避敌方威胁和环境障碍，安全地完成预定任务。随着无人机技术的飞速发展和广泛应用，在现代防空技术日臻完善的情况下，单纯地依靠手工操作无人机进行航路规划，已经无法有效地满足现代复杂任务的实际应用需求，于是一系列以解决最优化问题为核心的智能
优化算法脱颖而出，成为解决约束条件下无人机航路规划的重要方法。
[0003] 烟花算法由北大谭营教授于2010年在国际
群体智能大会上首次提出，其灵感来源于实际生活中烟花的爆炸过程，属于非
生物群体的智能优化方法，具有卓越的寻优性能，得到广大研究者的青睐，目前已广泛应用于聚类、网络规划、
图像识别、函数优化等领域问题的求解。烟花算法的基本思想在于：由多个个体(火花)组成群体，通过个体之间的相互协同使个体对环境的适应性变得越来越好，从而逼近最优解。在整个解空间中，烟花爆炸产生的火花会出现在以其为中心的邻域内，此过程即是对局部区域的一次搜索。在烟花算法中，为了平衡局部搜索能
力和全局探索能力，较好的烟花会在更小的爆炸范围内产生更多的火花，而较差的烟花会在更大的爆炸范围内产生少量的火花，于是较好的烟花负责在较优的区域进行更好的局部搜索，而较差的烟花负责在更大的爆炸范围内进行全局搜索，两者的协同搜索就是该算法所表现出的群体智能。
[0004]
粒子群优化算法是一种群体智能优化技术，1995年由Eberhart博士和Kennedy博士提出。该算法源于对
鸟群觅食行为的研究，是通过模拟鸟群觅食过程中迁徙和协同合作行为而提出的一种智能随机优化算法。其基本思想在于：各粒子从随机
位置出发，不断根据自身历史最优位置和最优粒子位置，调整搜索方向和速度，进而找到问题的最优解。粒子群优化算法以实现容易、
精度高、收敛快等优点引起了学术界的重视，并且在工程优化、
图像处理等诸多领域得到了广泛的应用。但是，粒子群优化算法在处理复杂问题时容易陷入局部最优解上而停止不动，即所谓的早熟收敛现象。
发明内容
[0005] 本发明提出了一种求解无人机约束航路规划的混合烟花粒子群协同算法(Hybrid FWPS Cooperation Algorithm,FWPSALC)。其目的在于提供一种有效的无人机航路规划方法，可以较快、较好地为无人机规划出一条远离威胁、路径最短的最优飞行路径，保障无人机以最短的时间、最安全的方式飞行到目标地点。本方法解决单一对象的路径规划问题，不涉及多对象之间的协同路径规划问题。
[0006] 混合算法中，采用两个种群并行独立搜索最优路径的方式。其中，一个种群采用改进的烟花算法进行搜索，而另一个种群则采用粒子群优化算法。当两个种群各自进行了一次搜索后，会进行信息的共享。烟花将通过爆炸在整个搜索空间中寻找到的最优解的大致区域提供给粒子，为粒子之后的搜索指引方向。而粒子则在接下来的
迭代过程中向该区域进行细致的局部搜索。如此，两个种群结合进行搜索进而得到规划问题的最优解。在整个搜索过程中，通过设置安全等级来划分安全路径和不安全路径，然后在一次次的迭代过程中，逐步提高安全等级，进而将最优路径的搜索范围限制在安全的路径中，以保证规划出的路径的安全性。
[0007] 具体技术方案如下：
[0008] 本发明提出的混合烟花粒子群协同算法兼顾了全局搜索能力和局部搜索能力，具有收敛速度快、鲁棒性强、寻优能力强等特点，能够有效地解决无人机约束航路规划问题，该方法的具体应用步骤如下：
[0009] 准备工作：无人机的起始点S、目标点T、威胁源个数H、各个威胁源的威胁等级及在飞行环境中的位置和威胁半径均为已知。算法的最大迭代次数为NCmax。
[0010] 步骤一：建立以ST为X'轴的旋转
坐标系X'-S-Y'，将起始点、目标点和障碍物的位置坐标转换到旋转坐标系中，计算规划空间的上下边界Smax、Smin，机器人规划出的路径为P＝{P0,P1,P2,…,PD,PD+1}，其中P0、PD+1分别表示起始点S和目标点T，Pi表示第i个路径点；
[0011] 步骤二：运用D条平行直线簇L1,L2,…,LD将ST垂直平分成D+1段，每段的距离为Δl＝||ST||/(D+1)，其中||ST||为ST的长度，路径点P1,P2,…,PD即对应位于直线L1,L2,…,LD上；
[0012] 步骤三：令迭代次数t＝1，初始化烟花，烟花个数为N，烟花维数为D，即随机生成N条无人机路径，每条路径中有D个路径点，烟花种群表示为向量组X(t)＝{X1(t),X2(t),...,XN(t)}，其中，Xi(t)＝{Xi1(t),Xi2(t),...,XiD(t)}表示第i个烟花，即对应第i条完整的路径，i＝1,2…,N，Xij(t)表示第i个烟花的第j个路径点，j＝1,2,...,D；
[0013] 步骤四：计算烟花种群中烟花对应路径P的路径长度fL(P)、威胁代价fT(P)、满意等级μ(P)，公式如下：
[0014]
[0015] 上式中||Pi+1-Pi||表示路径点Pi+1和Pi间的欧几里德距离；
[0016]
[0017] 上式中tj为第j个威胁源的威胁等级，P0.1,j,k为第j个威胁源对路径段PkPk+1上1/10点处的威胁概率，P0.3,j,k为第j个威胁源对路径段PkPk+1上3/10点处的威胁概率，依次类推。
[0018]
[0019] 上式中b为等级计算常数，其值为步骤四中所有烟花对应的路径的威胁代价的中位数。
[0020] 步骤五：初始化粒子群，设定粒子个数为Np，每个粒子对应一条完整的路径，每一个粒子的维数为D，即每条路径中有D个路径点，设定参数搜索范围、速度限制，随机确定每个粒子的初始位置和速度，即随机生成Np条无人机路径和对应的路径点变化量，并将粒子的位置存储在 中；
[0021] 步骤六：由公式(1)、(2)、(3)计算粒子群中各粒子对应路径P的路径长度fL(P)、威胁代价fT(P)、满意等级μ(P)；
[0022] 步骤七：将步骤四和步骤六中所有烟花和粒子中对应路径最短的个体保存在全局最优解xgbest(t)中，即目前搜索到的无人机最优路径，全局最优解xgbest(t)的初值即由初始化的烟花种群和初始化的粒子群计算得到；
[0023] 步骤八：更新安全等级α，公式如下：
[0024]
[0025] 上式中，ξ∈[0，100]为增长率参数,t为当前迭代次数；
[0026] 步骤九：计算烟花Xi(t)对应的爆炸火花数目Si和爆炸半径Ai，并生成爆炸火花Xj(t)＝[Xj1(t),Xj2(t),...,XjD(t)]，j＝1,2,...,Si，即对应无人机新的完整路径，Xjk(t)表示第j个爆炸火花的第k个路径点，k＝1,2,...,D；
[0027] 步骤十：将超出边界的爆炸火花映射到新位置，由公式(1)、(2)、(3)计算爆炸火花中各火花对应路径P的路径长度fL(P)、威胁代价fT(P)、满意等级μ(P)，所述的超出边界指火花上的路径点超出了规划空间的上下边界；
[0028] 步骤十一：通过等级比较方法选择烟花和爆炸火花中最优的N个个体，组成新的烟花群体V(t)＝{V1(t),V2(t),...,VN(t)}，所述的最优个体为通过等级比较方法得到的最优路径；
[0029] 步骤十二：对烟花Vi(t)进行变异操作，生成变异火花Vi′(t)，即对应无人机新的完整路径；
[0030] 步骤十三：将超出边界的变异火花映射到新位置，由公式(1)、(2)、(3)计算变异火花中各火花对应路径P的路径长度fL(P)、威胁代价fT(P)、满意等级μ(P)；
[0031] 步骤十四：通过等级比较方法选择烟花Vi(t)和变异火花Vi′(t)中更优的个体，i＝1，…，N，组成新的烟花群体U(t)＝{U1(t),U2(t),...,UN(t)}，其中最优烟花为Ubest(t)，最差烟花为Uworst(t)；
[0032] 步骤十五：更新粒子群中所有粒子的位置和速度，得到下一代粒子的位置和速度
[0033] 步骤十六：将超出边界的粒子映射到新位置，由公式(1)、(2)、(3)计算每个粒子对应路径P的路径长度fL(P)、威胁代价fT(P)、满意等级μ(P)，所述的超出边界的粒子指粒子上的路径点超出了规划空间的上下边界，或粒子上的路径点变化量超出了设置的速度范围；
[0034] 步骤十七：将更新后的第i个粒子与存储在xpbest(t)中的第i个粒子进行等级比较，选择更优的粒子保存在xpbest(t+1)中，i＝1,2,...,Np；
[0035] 步骤十八：通过等级比较方法选择下一代中的所有粒子x(t+1)和步骤14中的最差烟花Uworst(t)中最优的个体U′best(t)，并与步骤14的烟花群体U(t)中除Uworst(t)以外的其他烟花组成下一代烟花群体X(t+1)＝{X1(t+1),X2(t+1),...,XN(t+1)}；
[0036] 步骤十九：通过等级比较方法选择Ubest(t)、U′best(t)、xgbest(t)中最优的个体更新得到下一代全局最优解xgbest(t+1)；
[0037] 步骤二十：令t＝t+1，判断是否满足t>NCmax。如果不满足，返回步骤八进行迭代；如果满足，认为运算结束，输出全局最优解xgbest(t)。
[0038] 步骤二十一：将全局最优解xgbest(t)转换到X-O-Y坐标系中，输出规划的无人机最优路径。
[0039] 补充说明：
[0040] 步骤一中所述的坐标转换公式为：
[0041]
[0042] 上式中(xs,ys)为机器人起始位置，(x,y)、(x',y')分别为坐标系X-O-Y和X'-S-Y'中任意一个相同点的坐标位置，θ为X'轴与X轴的夹
角。
[0043] 步骤一中所述的规划空间上下边界为威胁源威胁区域中距离ST轴最远的点再向外侧延伸一定安全距离，其计算公式如下：
[0044]
[0045] 上式中ythreatj为坐标转换后第j个威胁源的Y'轴坐标，j＝1,2,...,H,Rj为第j个威胁源的威胁半径，max{}为取最大值函数，min{}为取最小值函数，Δd∈[5，50]为安全距离。
[0046] 计算路径威胁等级时，不同威胁源威胁概率计算如下：
[0047] (1)山峰威胁源威胁概率PP：
[0048]
[0049] 上式中d为无人机到威胁源的距离，RP为山峰威胁半径。
[0050] (2)雷达威胁源威胁概率PR：
[0051]
[0052] 上式中d为无人机到威胁源的距离，RRmax为雷达威胁半径。
[0053] (3)导弹威胁源威胁概率PM：
[0054]
[0055] 上式中d为无人机到威胁源的距离，RMmax为导弹威胁半径。
[0056] (4)高炮威胁源威胁概率PG：
[0057]
[0058] 上式中d为无人机到威胁源的距离，RGmax为高炮威胁半径。
[0059] 步骤五中，粒子群位置表示为向量组 其中，xi(t)＝[xi1(t),xi2(t),...,xiD(t)]表示第i个粒子的位置，即对应第i条完整路径的位置，i＝1,2…,Np，xij(t)表示第i个粒子第j个路径点的位置，j＝1,2,...,D，速度表示为向量组其中，vi(t)＝[vi1(t),vi2(t),...,viD(t)]表示第i个粒子的速度，
即第i条路径的变化量，vij(t)表示第i个粒子第j个路径点的变化量，j＝1,2,...,D，参数搜索范围为[Smin,Smax]，粒子速度最大值和最小值计算公式如下：
[0060] vmax＝0.01×(Smax-Smin) (11)
[0061] vmin＝-0.01×(Smax-Smin) (12)
[0062] 步骤九中，爆炸火花数目Si计算公式如下：
[0063]
[0064]
[0065] 公式(13)和(14)中，N1为对应路径安全的烟花个数，所述的安全路径指烟花Xi的满意等级μ(Xi)≥α的路径，N2为对应路径不安全的烟花个数，所述的不安全路径指烟花Xi的满意等级μ(Xi)<α的路径，β∈[0，1]为比例参数，M∈[10,100]是调整生成爆炸火花个数的常数，fLmax＝max(fL(Xi))为对应路径安全的烟花种群中最大路径长度值，i＝1,2,…,N1，μmin＝min(μ(Xi))为对应路径不安全的烟花种群中最小满意等级值，i＝1,2,…,N2，ε是机器最小量，为2-52，用于避免除零操作。
[0066] 当烟花Xi的满意等级μ(Xi)≥α时，即对应的路径安全时，其爆炸火花数目用公式(13)计算；当烟花Xi的满意等级μ(Xi)<α时，即对应的路径不安全时，其爆炸火花数目用公式(14)计算。为防止烟花Xi生成的爆炸火花过多或过少，对采用公式(13)或(14)计算出的爆炸火花数目Si，用如下公式进行限制：
[0067]
[0068] 上式中，round()是根据四舍五入原则取整的函数；a＝0.6、b＝0.05是两个常数。
[0069] 步骤九中，爆炸火花的爆炸半径Ai计算公式如下：
[0070]
[0071]
[0072] 公式(16)和(17)中，N1为对应路径安全的烟花个数，N2为对应路径不安全的烟花个数，是调整爆炸半径大小的常数，fLmin＝min(fL(Xi))为对应路径安全的烟花种群中最小路径长度值，i＝1,2,…,N1，μmax＝max(μ(Xi))为对应路径不安全的烟花种群中最大满意等-52级值，i＝1,2,…,N2，ε是机器最小量，为2 ，用于避免除零操作。
[0073] 当烟花Xi的满意等级μ(Xi)≥α时，即对应的路径安全时，其爆炸半径用公式(16)计算；当烟花Xi的满意等级μ(Xi)<α时，即对应的路径不安全时，其爆炸半径用公式(17)计算。为防止烟花Xi爆炸半径过小，对采用公式(16)或(17)计算出的爆炸半径Ai，用如下公式进行限制：
[0074]
[0075] 上式中，Amin∈[0,40]为烟花最小爆炸半径。
[0076] 烟花生成爆炸火花，其中，爆炸火花中第k个路径点的计算公式如下：
[0077] Xjk＝Xik+rand(-1,1)×Ai (19)
[0078] 上式中，rand()为取区间内随机值的函数。
[0079] 所述的火花映射规则如下：
[0080] Xjk＝Smin+|Xjk|％(Smax-Smin) (20)
[0081] 上式中，％为取模运算。
[0082] 所述的等级比较方法具体如下：
[0083] 迭代过程中，若两待比较个体R1、R2，对应路径长度f1、f2，路径满意等级μ1、μ2，安全等级α，采用等级比较的方法判断个体的优劣，具体方法如下：
[0084] (1)若μ1≥α、μ2≥α，且f1≤f2，选择R1为更优个体；
[0085] (2)当不满足(1)时，若μ1＝μ2，且f1≤f2，选择R1为更优个体；
[0086] (3)当不满足(1)和(2)时，若μ1>μ2，选择R1为更优个体；
[0087] 步骤十二中，变异火花生成公式如下：
[0088] Vi′＝Vi+rand(-1,1)×(Vb-Vr) (21)
[0089] 上式中，i＝1,2,...,N，Vb为烟花群体V(t)中最优的烟花，最优烟花即为最优路径，Vr为已生成的变异烟花和除Vi外V(t)中还未变异的烟花中的随机一个烟花；
[0090] 步骤十五中，粒子速度和位置更新公式如下所示：
[0091] vi(t+1)＝ω·vi(t)+c1·rand(0,1)·(xpbest,i(t)-xi(t))+c2·rand(0,1)·(xgbest(t)-xi(t)) (22)
[0092] xi(t+1)＝xi(t)+vi(t+1) (23)
[0093] 公式(22)中，i＝1,2,...,Np，c1∈[0,2]、c2∈[0,2]为
加速度常数，ω∈[0,1]为惯性权重，rand()为取区间内随机值的函数。
[0094] 步骤十六中，所述的粒子映射规则如下：
[0095]
[0096]
[0097] 上式中，i＝1,2…,Np，j＝1,2,...,D。
[0098] 有益效果：
[0099] 本发明提出了一种求解无人机约束航路规划的混合烟花粒子群协同算法。所提出的方法采用双种群并行搜索最优解的形式，一个种群采用改进的烟花算法，另一个种群采用粒子群算法，最后通过信息共享的方式找到无人机的最优航路。粒子群算法具有很好的局部搜索能力，弥补了烟花算法局部开发能力的不足，但是容易陷入局部最优解。而改进的烟花算法通过新的变异方式，得到了多样性较高的个体，提高了烟花群体的多样性，使得算法具备更好的全局搜索能力，进而有效地避免陷入局部最优解，弥补了粒子群算法的不足。双种群协同合作使算法具备了更优的全局搜索能力和局部搜索能力，同时采用等级比较的方法，较好的处理了约束条件，从而规划出更优的无人机飞行路径。混合烟花粒子群协同算法具有很好的优化性能，收敛精度高、鲁棒性强，且能规划出较短的安全路径，是一种有效的无人机航路规划方法。
附图说明
[0100] 图1是坐标转换示意图。
[0101] 图2是混合烟花粒子群协同算法
流程图。
[0102] 图3是混合烟花粒子群协同算法路径规划图。
[0103] 图4是混合烟花粒子群协同算法路径规划结果得统计数据表。
[0104] 图5是平均路径长度进化曲线图。
具体实施方式
[0105] 下面结合具体实施过程对本发明进行进一步描述，以图3的环境为例进行说明：
[0106] 准备工作：设无人机的飞行速度v＝1m/s。参数取值如下：最大迭代次数NCmax＝500；安全距离Δd＝20；增长率参数ξ＝15；爆炸火花粒子常数M＝30；比例参数β＝0.8；烟花最小爆炸半径Amin＝20；加速度常数c1＝0.8、c2＝0.8；惯性权重ω＝0.9。设起始点坐标为S＝[50,50]，目标点坐标为T＝[950,550]，烟花和粒子维度D＝25。
[0107] 步骤一：建立以ST为X'轴的旋转坐标系X'-S-Y'，将起始点、目标点和障碍物的位置坐标转换到旋转坐标系中，坐标转换公式如下：
[0108]
[0109] 上式中(xs,ys)为机器人起始位置，(x,y)、(x',y')分别为坐标系X-O-Y和X'-S-Y'中任意一个相同点的坐标位置，θ为X'轴与X轴的夹角。
[0110] 计算规划空间的上下边界Smax、Smin，公式如下：
[0111]
[0112] 上式中ythreatj为坐标转换后第j个威胁源的Y'轴坐标，j＝1,2,...,8，Rj为第j个威胁源的威胁半径，max{}为取最大值函数，min{}为取最小值函数。
[0113] 步骤二：运用25条平行直线簇L1,L2,…,L25将ST垂直平分成26段，每段的距离为Δl＝39.6；
[0114] 步骤三：令迭代次数t＝1，初始化烟花，烟花个数为8，烟花维数为25，烟花表示为向量组X(t)＝{X1(t),X2(t),…,X8(t)}，其中，Xi(t)＝[Xi1(t),Xi2(t),…,Xi25(t)]表示第i个烟花，i＝1,2…,8，Xij(t)表示第i个烟花的第j个路径点，j＝1,2,...,25；
[0115] 步骤四：计算烟花群体中烟花对应路径的路径长度、威胁代价、满意等级，公式如下：
[0116] 路径P的路径长度fL(P)计算公式如下：
[0117]
[0118] 上式中||Pi+1-Pi||表示路径点Pi+1和Pi间的欧几里德距离。
[0119] 路径P的威胁代价fT(P)计算公式如下：
[0120]
[0121] 上式中tj为第j个威胁源的威胁等级，P0.1,j,k为第j个威胁源对路径段PkPk+1上点处的威胁概率。
[0122] 不同威胁源威胁概率计算如下：
[0123] (1)山峰威胁源威胁概率PP：
[0124]
[0125] 上式中d为无人机到威胁源的距离，RP为山峰威胁半径。
[0126] (2)雷达威胁源威胁概率PR：
[0127]
[0128] 上式中d为无人机到威胁源的距离，RRmax为雷达威胁半径。
[0129] (3)导弹威胁源威胁概率PM：
[0130]
[0131] 上式中d为无人机到威胁源的距离，RMmax为导弹威胁半径。
[0132] (4)高炮威胁源威胁概率PG：
[0133]
[0134] 上式中d为无人机到威胁源的距离，RGmax为高炮威胁半径。
[0135] 路径P的满意等级μ(P)计算公式如下：
[0136]
[0137] 上式中b为等级计算常数，其值为步骤四中所有烟花对应的路径的威胁代价的中位数。
[0138] 步骤五：初始化粒子群，设定粒子个数为15，每个粒子对应一条完整的路径，每一个粒子的维数为25，粒子群位置表示为向量组x(t)＝[x1(t),x2(t),...,x15(t)]，其中，xi(t)＝[xi1(t),xi2(t),...,xi25(t)]表示第i个粒子的位置＝，i＝1,2…,15，xij(t)表示第i个粒子第j个路径点的位置，j＝1,2,...,25，速度表示为向量组v(t)＝[v1(t),v2(t),...,v15(t)]，其中，vi(t)＝[vi1(t),vi2(t),...,vi25(t)]表示第i个粒子的速度，vij(t)表示第i个粒子第j个路径点的变化量，j＝1,2,...,25，设定参数搜索范围[Smin,Smax]、速度限制vmin、vmax，计算公式如下：
[0139] vmax＝0.01×(Smax-Smin) (35)
[0140] vmin＝-0.01×(Smax-Smin) (36)
[0141] 随机确定每个粒子的 初始位置和速度，并将粒子的位置存储在中；
[0142] 步骤六：由公式(28)～(34)计算粒子群中各粒子对应路径的路径长度、威胁代价、满意等级；
[0143] 步骤七：将步骤四和步骤六中所有烟花和粒子中对应路径最短的个体保存在全局最优解xgbest(t)中，即目前搜索到的无人机最优路径；
[0144] 步骤八：更新安全等级α，更新公式如下所示：
[0145]
[0146] 步骤九：计算烟花Xi(t)对应的爆炸火花数目Si，计算公式如下：
[0147]
[0148]
[0149]
[0150] 公式(38)和(39)中，N1为对应路径安全的烟花个数，N2为对应路径不安全的烟花个数，fLmax＝max(fL(Xi))为对应路径安全的烟花种群中最大路径长度值，i＝1,2,…,N1，μmin＝min(μ(Xi))为对应路径不安全的烟花种群中最小满意等级值，i＝1,2,…,N2，ε是机器最小量，为2-52，用于避免除零操作。当烟花Xi的满意等级μ(Xi)≥α时，即对应的路径安全时，其爆炸火花数目用公式(38)计算；当烟花Xi的满意等级μ(Xi)<α时，即对应的路径不安全时，其爆炸火花数目用公式(39)计算。为防止烟花Xi生成的爆炸火花过多或过少，对采用公式(38)或(39)计算出的爆炸火花数目Si，用公式(40)进行限制。
[0151] 计算烟花Xi(t)对应的爆炸半径Ai，计算公式如下：
[0152]
[0153]
[0154]
[0155] 公式(41)和(42)中， 是调整爆炸半径大小的常数，N1为对应路径安全的烟花个数，N2为对应路径不安全的烟花个数，fLmin＝min(fL(Xi))为对应路径安全的烟花种群中最小路径长度值，i＝1,2,…,N1，μmax＝max(μ(Xi))为对应路径不安全的烟花种群中最大满意等级值，i＝1,2,…,N2，ε是机器最小量，为2-52，用于避免除零操作。当烟花Xi的满意等级μ(Xi)≥α时，即对应的路径安全时，其爆炸半径用公式(41)计算；当烟花Xi的满意等级μ(Xi)<α时，即对应的路径不安全时，其爆炸半径用公式(42)计算。为防止烟花Xi爆炸半径过小，对采用公式(41)或(42)计算出的爆炸半径Ai，用公式(43)进行限制。
[0156] 生成爆炸火花Xj(t)＝[Xj1(t),Xj2(t),...,XjD(t)]，j＝1,2,...,Si，Xjk(t)表示第j个爆炸火花的第k个路径点，k＝1,2,...,25，生成公式如下：
[0157] Xjk＝Xik+rand(-1,1)×Ai (44)
[0158] 步骤十：将超出边界的爆炸火花映射到新位置，映射公式如下所示：
[0159] Xjk＝Smin+|Xjk|％(Smax-Smin) (45)
[0160] 由公式(28)～(34)计算爆炸火花中各火花对应路径的路径长度、威胁代价、满意等级；
[0161] 步骤十一：通过等级比较方法选择烟花和爆炸火花中最优的8个个体，组成新的烟花群体V(t)＝{V1(t),V2(t),...,V8(t)}，采用等级比较的方法判断个体的优劣，具体方法如下：
[0162] 若两待比较个体R1、R2，对应路径长度f1、f2，路径满意等级μ1、μ2，安全等级α，[0163] (1)若μ1≥α、μ2≥α，且f1≤f2，选择R1为更优个体；
[0164] (2)当不满足(1)时，若μ1＝μ2，且f1≤f2，选择R1为更优个体；
[0165] (3)当不满足(1)和(2)时，若μ1>μ2，选择R1为更优个体；
[0166] 步骤十二：对烟花Vi(t)进行变异操作，生成变异火花Vi′(t)，变异公式如下所示：
[0167] Vi′＝Vi+rand(-1,1)×(Vb-Vr) (46)
[0168] 上式中，i＝1,2,...,8，Vb为烟花群体V(t)中最优的烟花，Vr为已生成的变异烟花和除Vi外V(t)中还未变异的烟花中的随机一个烟花。
[0169] 步骤十三：由公式(45)将超出边界的变异火花映射到新位置，由公式(28)～(34)计算变异火花中各火花对应路径的路径长度、威胁代价、满意等级；
[0170] 步骤十四：通过等级比较方法选择烟花Vi(t)和变异火花Vi′(t)中更优的个体，组成新的烟花群体U(t)＝{U1(t),U2(t),…,U8(t)}，其中最优烟花为Ubest(t)，最差烟花为Uworst(t)；
[0171] 步骤十五：更新粒子群中所有粒子的位置和速度，得到下一代粒子的位置和速度 更新公式如下所示：
[0172] vi(t+1)＝0.9·vi(t)+0.8·rand(0,1)·(xpbest,i(t)-xi(t))+0.8·rand(0,1)·(xgbest(t)-xi(t)) (22)
[0173] xi(t+1)＝xi(t)+vi(t+1) (23)
[0174] 公式(22)中，i＝1,2,...,15，rand()为取区间内随机值的函数。
[0175] 步骤十六：将超出边界的粒子映射到新位置，映射公式如下所示：
[0176]
[0177]
[0178] 上式中，i＝1,2…,15，j＝1,2,...,25。
[0179] 由公式(28)～(34)计算每个粒子对应路径的路径长度、威胁代价、满意等级；
[0180] 步骤十七：将更新后的第i个粒子与存储在xpbest(t)中的第i个粒子进行等级比较，选择更优的粒子保存在xpbest(t+1)中，i＝1,2,...,Np；
[0181] 步骤十八：通过等级比较方法选择下一代中的所有粒子x(t+1)和Uworst(t)中最优的个体U′best(t)，并与烟花群体U(t)中除Uworst(t)以外的其他烟花组成下一代烟花群体X(t+1)＝{X1(t+1),X2(t+1),…,XN(t+1)}；
[0182] 步骤十九：通过等级比较方法选择Ubest(t)、U′best(t)、xgbest(t)中最优的个体更新得到下一代全局最优解xgbest(t+1)；
[0183] 步骤二十：令t＝t+1，判断是否满足t>NCmax。如果不满足，返回步骤八进行迭代；如果满足，认为运算结束，输出全局最优解xgbest(t)。
[0184] 步骤二十一：将全局最优解xgbest(t)转换到X-O-Y坐标系中，输出规划的无人机最优路径。
[0185] 图3显示了通过上述计算步骤得到的无人机航路规划结果，图4表示了混合烟花粒子群协同算法50次独立航路规划的数值结果，图5为混合烟花粒子群协同算法50次独立航路规划的平均路径长度进化曲线。由图3可以看出采用混合烟花粒子群协同算法规划出的无人机航路既安全又平滑，是一条很好的无人机飞行路径。图4中第二列为50次规划中规划出的最
短路径的长度，第三列为50次规划中所有路径长度的平均数值，第四列为50次规划中所有路径长度的中位数，第五列为50次规划中路径长度的标准差，第六列为50次规划中规划出的最长路径的长度，第七列为50次规划的成功率，规划成功即无人机能安全地从起始点到目标点，第八列为50次规划中算法规划一次的平均时间。由图4中的数据可以得出：采用混合烟花粒子群协同算法进行无人机航路规划，可以在允许的规划时间里，规划出较优的无人机飞行路径，且算法鲁棒性最强，成功率高。因此本发明所提出的混合烟花粒子群协同算法是一种有效的无人机航路规划方法，可以稳定地规划出较优的无人机飞行路径，并保障无人机可以安全运行到目标点。
