基于自适应黏菌算法优化的无人机三维路径规划

黄 鹤, 高永博, 茹 锋, 杨 澜, 王会峰

(长安大学 a. 电子与控制工程学院; b. 西安市智慧高速公路信息融合与控制重点实验室;c. 信息工程学院,西安 710064)

无人机具有生存能力强和机动性能优越的特点,被应用于民生和国防等诸多领域[1-2].设计高效合理的无人机路径规划方案,可以保证无人机在执行任务时有效躲避各种威胁区域,顺利到达目的地.近年来许多国内外学者对无人机的路径规划方法进行了大量的研究,对此提出了一些路径规划算法,主要可以分为两大类:① 传统算法,如人工势场法[3]和A*算法[4]等;② 智能算法,如粒子群算法[5]、灰狼算法[6]和蝠鲼觅食算法[7]等.采用群体智能优化算法解决无人机路径规划的问题是目前的研究热点,相比于传统算法,智能算法具有寻优速度快、稳定性强的优点.在此基础上,王翼虎等[8]利用改进粒子群优化算法解决无人机路径规划的问题,生成一条航程较短的飞行路径,收敛速度更快,但收敛精确度不高.黄书召等[9]为无人机成功通过威胁区域提出了一种改进遗传算法的无人机三维路径规划,提升了搜索性能,但威胁区域考虑得比较单一,实际应用效果并不理想.吴坤等[10]利用贪婪优化策略改进鲸鱼算法,通过扩大搜索路径提升全局寻优能力,收敛速度较慢.上述研究算法虽然实现了无人机的路径规划,但算法的收敛速度和寻优精度有限.2020年,Li等[11]根据黏菌在觅食过程探索路径的行为提出了黏菌算法(Slime Mould Algorithm,SMA),特点是寻优速度快、算法简单、模型易修改,适合应用在无人机路径规划中,但SMA自身存在寻优精度不高和容易陷入局部最优的缺点,仍需进一步改进.因此,本文提出了一种基于改进的自适应黏菌算法(GSMA),保证了寻优速度和寻优精度,又避免陷入局部最优,应用在无人机中能够提高其路径搜索能力.

1 无人机三维路径规划建模

在无人机三维路径规划的过程中,根据任务、威胁源情况和地形等实际因素,构建了三维地形和威胁源约束.将这些威胁源结合自身约束条件建模,等效出三维环境下总的地形约束.

1.1 地形约束

地形环境是无人机路径规划首先要考虑的问题,不同的地貌环境对无人机的威胁不同.三维地形可大致分为山地、丘陵和平原,其中山地条件中的山峰对无人机的正常飞行影响最大,当然也不能忽略海拔的影响.山峰建模坐标Z(x,y)如下所示:

(1)

式中:坐标(x,y)和(x0,y0)分别表示三维地形下山峰在地平面上和中心点的位置;h表示山峰的高度;λ1、λ2表示山峰的倾斜度.无人机路径规划过程中的地形威胁代价函数表示如下:

(2)

fHj=hj

(3)

要充分保证无人机的飞行安全,在路径规划中除了要考虑地形约束的影响,还需要考虑飞行的边界和最高距离.实验中,假设无人机的飞行范围表示为(Xmin,Ymin)=(0, 0) km和(Xmax,Ymax)=(100, 100) km,探索路径过程中的最高飞行高度为Zmax=5 km.

1.2 威胁模型约束

无人机在飞行过程中的威胁主要来自雷达、电磁威胁和导弹截击等,建模如下.

(1) 雷达威胁.

雷达主要通过电磁波探测飞行物的距离和速度信息,对无人机的安全通过威胁很大.雷达侦测飞行物的方程可以简化为圆锥形探测模型,描述如下:

P0=(XM-X0)2+(YM-Y0)2=

(4)

式中:(X0,Y0,Z0)表示雷达中心地面坐标;(XM,YM,ZM)表示雷达威胁区最大位置坐标;h0代表雷达威胁源的高度;r0为圆锥模型底面半径.

(2) 电磁威胁.

电磁威胁是无人机路径规划的又一主要威胁.一般情况下,可以用半球形等价电磁威胁区域.电磁威胁函数模型PC可表示为

(5)

式中:X、Y、Z分别指三维平面的X轴、Y轴和Z轴;R代表电磁威胁半球形的区域半径;α和β分别为三维平面上Z轴正向与电磁威胁半径的夹角以及威胁半径在水平面上与X轴正向夹角.

(3) 导弹威胁.

敌方导弹威胁到无人机的正常工作,飞行过程中必须躲避.导弹威胁区域Pd1、攻击距离Pd2和威胁总代价PD模型表示如下:

(6)

(7)

PD=Pd1+Pd2

(8)

式中:kd和kg分别表示导弹击中无人机的概率和导弹的加速度参数;rd代表航迹点到导弹威胁中心的距离;rc表示导弹影响范围半径;rg表示无人机和目标之间的直线距离.

(4) 禁飞区威胁.

无人机在飞行过程中还会受到恶劣天气和复杂环境的影响,被称为禁飞区.禁飞区的代价表示为

PF=KF

(9)

式中:KF为禁飞区的威胁代价.

将威胁源等效为地形模型,假设距离威胁源中心越近,威胁代价越大,地形越高,反之则越低.因此,可以将威胁源等效为

Qjp=

(10)

式中:Kthr为威胁源修正系数;Rmax,p表示第p个威胁源的最大半径;第p个威胁源中心的水平坐标为(xp,yp).禁飞区的威胁区域的代价fjp表示如下:

(11)

式中:rj,p为航迹点j到威胁源中心p的距离.

无人机在进行执行任务的同时与地形威胁山源保持安全距离,避免产生碰撞,同时规避雷达、电磁和导弹威胁,保障无人机飞行安全.

1.3 无人机自身约束

无人机抽象为一个质点模型,无需考虑自身的质量和形状.在路径规划的过程中,无人机不仅受到外部条件威胁,也受限于自身约束条件,如自身转弯角度α0、向上爬升角度β0和燃油消耗(用飞行航程表示).设最大转弯角度αmax、向上最大爬升角度βmax,这些自身物理约束条件可以表示为

(12)

(13)

(14)

式中:Qh和Qv分别为无人机的转弯角和爬升角系数;Jh和Jv分别为转弯角度和向上爬升角度的代价函数;JLj为航迹点j对应的航程;lj为无人机的飞行航程.综合上述所有的代价函数,无人机航迹点总约束代价表示为

fJj=Jh+Jv+JLj

(15)

1.4 航迹代价函数

将上述地形约束、高程代价、威胁模型约束及无人机自身物理约束的代价加权综合起来,得到无人机总的代价函数,用F(Rj)表示为

F(Rj)=

(16)

式中:d为航迹点数目;σ威胁源总数;ω1、ω2、ω3、ω4为各约束代价的权重.

2 改进的黏菌觅食算法

2.1 黏菌算法

黏菌是一种真核生物,SMA模拟了黏菌在捕食过程中的行为变化,利用黏菌的觅食过程进行路径规划.

2.1.1黏菌觅食位置的更新 黏菌在觅食过程中,空气中食物的浓度越大,黏菌自身细胞质流动越快,进而觅食成功.黏菌觅食的过程可用函数表达式表示,位置更新公式如下:

Xnew=r(BU-BL)+BL,r

(17)

Xnew=

(18)

(19)

p=tanh(S(i)-DF)

(20)

式(17)和(18)表示黏菌觅食的过程,BU和BL规定了黏菌搜索食物范围的上下边界;c为常量,文中取0.03;r的取值范围是在0至1区间内的随机值;Vb参数在-a到a区间内取值;Vc参数在觅食过程中线性减少至0;t表示黏菌的迭代次数;Xb代表食物浓度最高的位置;X为黏菌的位置;XA和XB代表黏菌的两个任意位置;W为黏菌自身的质量;Tmax表示最大迭代次数.黏菌位置X可以随着其获得的最优位置的变化而更新,同时黏菌位置也可以根据Vb、Vc和W参数的微调而更新[12].式(20)中的S(i)为黏菌的适应度,i表示种群的位置,而DF为黏菌的最佳适应度.此外,黏菌质量参数W在不同食物浓度的函数表达式表示如下:

(21)

Sr=sort(S)

(22)

式中:N/2表示前一半黏菌的数量;bF和wF分别为在迭代过程中的最优适应度和最差适应度;sort表示对种群的适应度进行升序排列;Sr表示呈递增型的适应度序列.

2.1.2算法流程 SMA流程如下.

步骤1初始化种群,设定参数.

步骤2计算种群适应度值并排序.

步骤3利用式(17)和(18)更新黏菌种群位置.

步骤4计算种群适应度值,并更新最优位置,计算当前最优位置.

步骤5判断是否满足最优条件,若满足则输出最优结果,否则执行步骤2～5直至满足最优.

2.2 基于改进的SMA的无人机三维路径规划

2.2.1航迹编码 无人机飞行过程中的任一条路径被随机的一只黏菌所定义,而飞行路径由许多个航迹点连接而成的线组成,并且每个航迹点均具有三维空间属性(x,y,z).路径规划过程中,要确定航迹点数,然后将z方向的坐标固定不变,即将三维等效为二维问题解决,则只需要对x和y方向坐标进行寻优.

2.2.2改进的Logistic映射 由于SMA存在过度依赖种群初始位置的问题,以随机方式初始种群易使得种群分布不均匀,影响算法的求解精度.混沌运动具有遍历性和不重复性等特点,适合用来生成种群的初始位置[12].大多群体智能优化算法初始种群映射方式使用的是Tent映射和Logistic映射,Tent映射在0～2的取值范围内存在着小周期现象、容易陷入不动点的问题,而传统Logistic映射对初始值条件要求严格,使得初始化黏菌种群的分布不均且多样性不足.据此,设计一种改进的Logistic混沌映射,以增加种群的多样性和均匀分布.具体表达如下:

xtemp=mod(a/xn+λxn(1-xn),b)

(23)

(24)

式中:mod为取模函数;xtemp为中间变量;mod(a/xn,b)的取值范围为[0, 1);λ为常数,文中取0.3;b为(0, 1)中的随机数,文中设定b=0.7;最后输出序列xn的范围是(0, 1].式(23)将不会出现不定值,且在a不同取值下仍具有较好的随机均匀性.

2.2.3非线性自适应权重因子 权重因子是SMA的核心参数,设计自适应权重因子可以用于改善算法的搜索能力[13]. SMA中参数是线性变化的,搜索能力较差且搜素范围局受到限制,寻优能力较弱.为解决这个问题,设计了一种非线性自适应权重因子w,黏菌位置更新公式如下:

Xnew=

(25)

(26)

式中:wmax和wmin分别为惯性权重因子的最大值和最小值,本文经过大量实验取值为wmax=0.82,wmin=0.01.自适应权重因子w系数变化如图1所示.

图1 w系数变化图

在迭代前期的权重因子w快速增大,扩大种群的搜索范围,w在迭代后期增速放缓,有利于黏菌觅食过程中跳出局部最优,提高了算法的精度,有益于寻找到更合适的位置[14].自适应权重变化是指黏菌觅食过程中,黏菌位置会根据黏菌的迭代次数的变化而更新.GSMA设计的非线性自适应惯性权重因子,增强了搜索能力,同时避免了陷入局部最优.

2.2.4自适应柯西变异策略 由于地形情况较为复杂,在黏菌觅食过程中,对最优位置进行柯西变异[15],可以提高迭代后期的搜素能力.柯西分布的特点是中间的峰值比较小而两边范围比较广,这样可以使黏菌在当前范围产生更大的扰动.针对路径规划过程中出现陷入局部最优和搜索能力较弱的情况,标准柯西变异处理后优化效果有限,因此提出了一种自适应柯西变异策略,公式如下:

(27)

(28)

式中:θ为常数,在文中取10;rmax为黏菌个体之间的最大距离;F(γ;t)为柯西变异累计分布函数.自适应柯西变异在黏菌群体迭代过程中可以获得较多的最优值,保证平滑收敛至最优,增强GSMA的全局搜素能力,减少了黏菌最优值的动荡并提高了其寻优速度.自适应柯西变异下的黏菌位置更新如下所示:

Xbnew=Xnew+F(γ;t)

(29)

式中:Xnew为上一代的黏菌最优位置;Xbnew为经过自适应柯西变异更新后黏菌新的位置.

图2 GSMA收敛过程

当解位于Dmax区域时,算法的解个体会向着最优值xmax靠近,即方程的解x(t;x0,t0) 就会向着xmax靠近,此时解的过程表达式为

δ≤min{f(xmax)-f(x1),

f(xmax)-f(x2)}≤ε

(30)

(31)

因此,圆S(δ)的半径δ与初始状态t0无关,所以GSMA的平衡状态xmax是一致稳定的,即

(32)

综上所述,GSMA可以收敛到全局最优值.

2.2.6GSMA流程 GSMA流程如图3所示,具体流程如下.

图3 算法流程图

步骤1设置种群数量、最大迭代次数和相关参数.

步骤2利用改进Logistic混沌映射初始化种群,计算种群的最优和最差适应度并排序.

步骤3判断随机数r

步骤4计算适应度值,更新最优位置.

步骤5根据每次适应度的不同进行排序.

步骤6对当前位置进行自适应柯西变异,计算适应度值,更新最优位置.

步骤7是否到达最大迭代数,满足则输出最优值,不满足则重复上述步骤直至满足条件.

3 GSMA对比实验

3.1 测试函数

利用不同的测试函数对算法性能优劣进行比较,在不同复杂度的情况下验证GSMA寻优能力的强弱.测试函数具体描述如下所示.

(1) Rosenbrock函数.

(33)

式中:Rosenbrock函数搜索范围为[-30, 30],维度为30,全局最优值为0.

(2) Schwefel函数.

(34)

式中:Schwefel函数搜索范围为[-500, 500],维度为30,全局最优值为 12 569.5.

(3) Foxholes函数.

(35)

式中:Foxholes函数搜索范围为[-65.536, 65.536],维度为2,全局最优值为0.998.

(4) Kowalik函数.

(36)

式中:Kowalik函数搜索范围为[-5, 5],维度为4,全局最优值为0.000 3.

(5) Hartman 6函数.

(37)

式中:Hartman 6函数搜索范围为[0, 1],维度为6,全局最优值为 -3.332.

(6) Shekel 10函数.

(38)

式中:Shekel 10函数搜索范围为[-5, 5],维度为4,全局最优值为 -10.536 3.

3.2 实验验证

利用Rosenbrock、Schwefel、Foxholes、Kowalik、Hartman 6和Shekel 10共6种测试函数评估SMA、灰狼优化(Grey Wolf Optimizer, GWO)算法[16]、海鸥算法[17](Seagull Optimization Algorithm, SOA)以及GSMA的性能.经过30次测试后,适应度F随着迭代次数t变化如图4所示.适应度的最优值Fbest及均值Favg如表1所示.

表1 4种算法在测试函数上的实验对比

图4 各测试函数上的适应度变化曲线

由图4和表1可知,针对Rosenbrock函数,SMA收敛速度弱于GWO、SOA和GSMA,迭代80次时收敛并陷入局部最优;而另外3种算法均寻得最优,其中GSMA的收敛速度和寻优精度均最优.针对Schwefel函数,SMA和GWO过早收敛,SOA的寻优速度不如GSMA.经过改进后的GSMA寻优能力更好,优于GWO和SMA.针对Foxholes函数,GSMA迭代不足10次已经收敛,而GWO和SOA则陷入了局部最优.针对Kowalik函数,GWO和SOA迭代10次已收敛并陷入局部最优,SMA迭代90次时陷入局部最优,而GSMA迭代60次寻得最优,且收敛速度和寻优精度均为最优.针对Hartman 6函数,GWO、SOA和SMA分别迭代60、60和80次寻得最优,SMA的寻优精度更高,GSMA收敛速度则更快且寻得最优.针对Shekel 10函数,SOA过早陷入局部最优,SMA寻优精度和速度较低,而GSMA和GWO寻优精度大致相同,但收敛速度明显高于GWO和SMA.综上所述,SMA改进后提升了寻优速度和精度,更易跳出局部最优.

4 无人机突防仿真测试及分析

实验平台为Window10系统的计算机,CPU为AMD R5-5600H,频率为3.30 GHz、内存为16 GB,软件采用MATLAB 2021a.测试区域为100 km×100 km,高度为5 km,采用两种不同的仿真地形.在地形1中,无人机的起始与终止坐标分别为(1, 93, 0.2) km、(98, 8, 0.5) km,3座不同类型的雷达位置坐标分别为(32, 40, 0) km、(70, 58, 0) km、(40, 12, 0) km,侦测半径设定为9、11、14 km,电磁和导弹威胁半径分别为12、18 km.地形2的区域大小也为100 km×100 km,无人机的起始点和终点为(2, 3, 0.2) km、(99, 87, 0.5) km,3座不同类型的雷达位置坐标分别为(30, 40, 0) km、(69, 70, 0) km、(40, 11, 0) km,侦测半径设定为6、9、10 km,电磁和导弹威胁半径分别为10、16 km,飞行过程中的最大转弯角为60°.为提高无人机的生存和路径规划能力,两种地形山峰位置有所不同.仿真种群数量均为100,最大迭代次数也是100.建模的主要参数如表2所示,表中:v为飞行速度.两种地形中,无人机的三维路径的仿真结果如图5和图6所示,各算法仿真结果如表3和表4所示.

表2 主要建模参数

表3 地形1仿真结果统计

图6 地形2仿真结果图

地形中黄色波为不同类型的雷达探测区域,蓝色波表示电磁和导弹威胁区域.为对于仿真地形1,由图5(a)和5(b)可以看出,在SOA结果中,无人机

成功穿过了第1座山峰,但在第2座山峰处撞山且被敌方雷达侦测到,不能成功到达目的地;而GWO和SMA结果中,无人机虽然成功规避了山峰和雷达侦察,但没有绕开电磁威胁区域和导弹攻击区域,且路线过长,导致撞击山脉的概率增加,影响实用价值;相比之下,GSMA不仅成功规避了山脉和雷达侦测,而且选择了最优路线以更小的代价穿过地形.从表3可以看出,SMA、GWO、SOA以及GSMA的航程分别为166.23、159.38、162.17、133.54 km,最优适应度分别为30.86、34.09、36.72、29.43,GSMA的规划路径最短,转弯及爬升角代价最低,整体保持低空平稳飞行,航迹代价最低,优势明显.图5(c)的仿真结果表明,SOA和GWO在迭代次数为80时仍未寻得最优而陷入了局部最优,寻优性能相比于SMA较差,并不适合无人机的径规划;SMA的精度相对较低,迭代至90次陷入局部最优,且路径规划距离过长.

相比于地形1,图6中地形2的山峰更为密集,导弹和电磁覆盖威胁范围更广.从图6(a)和6(b)可以看出, SOA、GWO和SMA都规避了密集连片山脉,但经过了导弹威胁区域,降低了飞机的安全性且规划的路线过于漫长,实用性不高.GWO与山峰发生碰撞,未能安全到达目的地.而GSMA规避了雷达的侦测,选择最佳路线快速穿过地形2.相比之下,GSMA寻优精度和寻优速度都优于GWO及SMA,在迭代至80次时已经寻得最优.而此时GWO及SMA仍未收敛,SOA则迭代20次即陷入了局部最优.综上,GSMA能够更为有效地避开各种威胁地形,路径规划中的总代价也最小.

5 结语

提出了一种基于改进SMA的算法来解决无人机三维路径规划的问题.首先,建立相关的地形模型、各种威胁源和无人机约束模型,并建立代价函数.其次,采用混沌映射对种群初始化,用以增加种群的多样性和随机分布,扩大了其搜索范围更易获得最优;在SMA中引入了非线性自适应权重因子,加快算法的收敛速度,很好地弥补了SMA自身存在的缺陷;加入柯西变异,使得种群在前期跳出局部最优,提高收敛速度.实验结果表明,GSMA在无人机探索路径规划过程中,能够在最短航程的条件下快速穿过危险区域和规避障碍物,应用价值明显.