目前，无人飞行器(unmanned aerial vehicle，UAV)越来越广泛地运用于军事和民用领域．作为UAV自主导航系统的核心，航迹规划是当前的研究热点．其目的是在给定的条件下，为UAV设计出从起始点到目标点的最优或满意的飞行航迹[1-2]．UAV航迹规划在本质上是一个复杂带约束的非确定性多项式难题(non-deterministic polynomial-hard，NP-hard)[1]．与传统的确定性算法相比，粒子群、蚁群和头脑风暴优化(brain storm optimization，BSO)算法等随机性算法更适合于求解复杂的NP-hard优化问题[2-3]．传统的群智能算法的群体智能性较低，存在着早熟收敛等诸多问题[2]．BSO算法作为一种较新的群智能算法，在解决复杂的优化问题方面具有一定的优势，被广泛地运用于解决车间调度等各种复杂优化问题[3]．BSO算法自提出以来，众多学者对其进行了研究并提出了多种改进的BSO算法．改进策略主要围绕种群的初始化方法、聚类方法、更新机制，以及与其他算法结合等方面[3-7]．然而，传统BSO及其相关改进算法仍存在容易陷入局部极值等问题．如何提升BSO算法的优化性能仍须进一步研究．本研究针对复杂环境中的UAV三维航迹规划问题，提出了基于改进BSO(quantum-behaved BSO，QBSO)[5-7]算法的航迹规划方法．与量子QBSO等算法相比，所提算法在解决UAV三维航迹规划问题方面具有更好的全局优化能力和更强的稳定性．1 UAV三维航迹规划数学模型1.1　环境建模和航迹编码如图1所示，战场环境中有待要攻击的目标和雷达威胁等威胁区以及敌方防空阵地等禁飞区．用三维圆柱描述威胁区和禁飞区．为便于计算，对规划地图进行了坐标旋转，并且只选取规划空间中一定大小的区域进行航迹规划．设CS'和CG'分别是出发点CS和目标点CG在原平面xoy的投影点．新坐标系以CS'为原点，以向量CS'CG'⃗所在的直线为x'轴，以x'轴在水平平面中的垂线为y'轴，以原坐标系的z轴为z'轴．10.13245/j.hust.240416.F001图1坐标轴旋转及航迹点的x'坐标等间隔采样示意图一条航迹X可以表示为三维空间中的一系列的航迹点X={CS,C1,C2,…,CD,CG}，(1)式中：Ci(i=1,2,...,D)为航迹中除出发点和目标点之外的航迹点，具有三维坐标(xi,yi,zi)；D为除出发点和目标点之外的航迹点的数量．如图1所示，为降低计算量，航迹中所有节点的x'轴坐标是对向量CS'CG'⃗等间隔采样得到的．在新坐标系中，第i个节点x'轴坐标为xi=iCS'CG'⃗ /(D+1)，其中⋅为向量的欧氏距离．可以看出：通过固定各节点的x'轴坐标的值，只须对D个航迹节点的y'轴和z轴坐标进行优化，从而将航迹规划问题的维数转化为2D．1.2　代价模型和约束条件UAV航迹的性能与航迹总长度等多种指标有关，同时UAV航迹规划受约束条件限制．这些约束条件可分为三类[1-2,8-10]：环境约束，主要包括敌方防空阵地等禁飞区约束，可将禁飞区用威胁系数为无穷大的威胁区描述；最大拐弯角约束等UAV自身机动性能约束；出发角度约束等任务约束．假定UAV以GPS导航方式执行作战任务，其作战任务没有出发及到达角度约束．主要考虑禁飞区约束、UVA的最大转弯角约束、最大爬升/下滑角约束、最大及最小离地约束等约束条件．这些约束条件已经涵盖了本研究设定任务场景下UAV航迹规划主要的约束条件．航迹规划的目的就是在满足各种约束条件下，在规划空间内规划出一条总代价最小的航迹．将航迹代价定义为J(X)=∑i=14wiJi(X)，(2)式中：J1(X)~J4(X)为航迹的长度、威胁、机动和高度代价；w1~w4为各部分代价的权系数．航迹X的长度代价J1(X)是指整条航迹的总长度．设节点C0和CD+1分别代表出发点CS和目标点CG．J1(X)可表示为J(1X)=∑i=0DCiCi+1⃗．(3)UAV在飞行过程中应尽量远离防空导弹等敌方威胁．与文献[8-9]类似，航迹X的威胁代价与其穿越威胁区的长度有关．设Ci'和Ci+1'分别是航迹点Ci和Ci+1在水平平面的投影，Q为威胁区的个数．航迹段CiCi+1⃗的威胁代价T(Ci,Ci+1)定义为T(Ci,Ci+1)=∑k=1QρkLi,i+1k，(4)式中：Li,i+1k为Ci'Ci+1'⃗包含在第k个威胁区在水平平面内的投影长度；ρk为第k个威胁区的威胁系数．因此，威胁代价J2(X)可定义为X的所有航迹段的威胁代价之和J2(X)=∑i=0DT(Ci,Ci+1)．(5)J3(X)定义为所有节点Ci(i=0,1,…,D)处的转弯代价Fi和爬升/下滑代价Si之和J3(X)=∑i=0D(Fi+Si)，(6)式中：Fi=αi  (αi≤αmax),E      (αiαmax);Si=φi  (φi≤βmax),E   (φiφmax); (7)αmax和βmax为UAV的最大转弯角和最大爬升/下滑角；E为惩罚因子；αi和φi为航迹点Ci处的转弯角和爬升/下滑角．在地形跟随过程中，UAV与地面之间必须保证最小安全间隙hmin，不能超过最大离地高度hmax．设zi为节点Ci=(xi,yi,zi)的高度，Hi为(xi,yi)处的地形高度．J4(X)可表示为：J4(X)=∑i=1Ddi；(8)di=zi-Hi-hmin    (hmin≤zi-Hi≤hmax),E   (zi-Hihmax,   zi-Hihmin). (9)2 改进QBSO算法三维航迹规划2.1　改进QBSO算法基于改进QBSO算法的航迹规划流程如图2所示，Tx为协同进化代数阈值；G为最大进化代数．为更好地平衡全局搜索和局部搜索，采用协同进化机制且整个种群的大小在进化前期和后期不同．在进化前期，两个种群独立进化且种群之间没有信息交流，从而可保证种群之间的差异性并提高算法全局搜索能力．在进化后期，每个种群中较优的一半个体合并为一个新种群．新种群按照QBSO算法的进化机制继续进行进化，从而提高算法的收敛速度和收敛精度．为进一步提升算法的优化能力，提出了一种待变异个体产生方式．10.13245/j.hust.240416.F002图2改进QBSO算法的航迹规划流程2.2　个体的编码及种群的初始化种群中每个个体Xi的维数为2D，可表示为一个维数为2D的向量．Xi=[xi1,xi2,...,xi2D]，式中，Xi的前1/2维和后1/2维分别代表航迹Xi上D个节点在新坐标系下的纵坐标和高度．设Uij为航迹Xi上第j个节点在的空间坐标，Uij=jCS'CG'⃗ /(D+1),xij,xij+D(j=1,2,...,D),根据式(2)～(9)可计算航迹Xi的总代价．对个体Xi进行初始化时，Xi上每个航迹点在新坐标系下的纵坐标和高度按如下方式随机生成xij=(r1-0.5)W；xij+D=Hij+hmin+r2(hmax-hmin)，式中：r1和r2为区间[0,1]内的随机数；W为规划区域的宽度；Hij为点jCS'CG'⃗ /(D+1),xij处的地形高度．为提高算法的全局搜索能力，按照上述初始化方法初始化两个同等大小的子种群．2.3　子种群独立进化在进化前期，每个子种群按照QBSO算法的原理独立进行进化，从而提高整个种群的多样性．在每个子种群的每一轮迭代过程中，首先将种群聚类为K个类别，大部分BSO算法的相关改进算法采用K-means算法对种群进行聚类[3-5,7,11]．近年来，也有研究者对种群聚类方法进行了研究，提出了多种聚类方法[6,12-14]，但各有优缺点．本研究采用K-means算法进行聚类(也可采用其他聚类方法)．此外，现有文献大都根据经验值确定K-means算法的类别数K[3-5,7,11]．文献[14]提出了一种用于确定最大适应度值聚类法类别数的自适应控制法．也可将K的取值问题归结为一个优化问题，通过运用智能优化算法对聚类数的指标函数进行优化，从而得到最优聚类个数[15-16]．接着，以一定概率pc产生随机个体取代某一个类的类中心，并以一定的方式产生待变异个体及新个体．在传统BSO算法的待变异个体产生方式的基础上，提出一种改进的待变异个体产生方式，即待变异个体还可以是3个类中心融合得到或者是3个类中的随机个体融合得到．所提改进的待变异个体产生方式增加了个体与其他类的交流，可以丰富种群的多样性，从而提升算法的全局搜索能力．如图3所示，每一子种群在每一代进化过程中，以下步骤产生待变异个体及新个体，p1~p4为各待变异个体产生方式的概率，r1~r5为区间[0,1]内的随机数．10.13245/j.hust.240416.F003图3改进的待变异个体产生方式的示意图当基于3个类产生待变异个体时，待变异个体的生成方式为Xiold=r1PXiold1+r2PXiold2+r3PXiold3，式中：Xiold为待变异个体；P=r1+r2+r3；Xiold1，Xiold2和Xiold3为分别来源于3个不同的类的个体．若当前进化代数小于阈值Tx，则重复上述迭代过程，否则按照下述的步骤进行后续进化．两个子种群独立进化Tx代后，首先对每个子种群中的个体进行排序，然后分别选取每个子种群中较优的(排名前50%)个体，并将这些较优的个体合并为一个新种群，最后新种群继续按照QBSO算法的进化机制和所提改进的待变异个体产生方式进行进化，直至达到算法终止条件，输出最优解．2.4　复杂度分析改进QBSO算法的时间复杂度主要由以下几部分组成[11]：a．初始化两个大小为N的种群，并计算每个个体的代价值的时间复杂度O(ND)；b．利用K-means聚类算法将种群聚为K类的时间复杂度为O(KND)，整个进化过程中聚类操作的时间复杂度为O(MKD)，M为代价函数计算总次数；c．在整个进化过程中，产生新个体并计算其代价值的时间复杂度O(MD)；d．当进化代数等于Tx时，对每个子种群进行排序的时间复杂度为O(Nlog2N)，则改进QBSO算法的时间复杂度为O(MKD)．BSO，QBSO和IBSO算法的时间复杂度也是O(MKD)．GBSO算法采用基于适应度值的分类方法，整个进化过程中其分类操作的时间复杂度为O(Mlog2N)．因此GBSO算法的时间复杂度为O(MD)(假设log2ND)．3 仿真实验分析选取战场环境比较有代表性的三组实验，对基于改进QBSO(Improved QBSO，IQBSO)算法与BSO、QBSO、全局最优BSO(Global-Best BSO，GBSO)[6]、改进BSO(Improved BSO，IBSO)[7]算法进行比较分析．实验的运行环境为Win10操作系统，CPU为intel 2.80 GHz，16 GiB RAM，编程环境为Matlab 2021a．实验使用了分辨率为373×373像素、大小为5 000×5 000像素的数字地形高程图和模拟生成的威胁数据．参照文献[3,5-8]中参数的设置方法，5种算法的共有参数的值均相同．综合考虑规划时间与规划精度并参照文献[8-10]，经过反复试验，参数设置如下：M=20 000，N=40，D=6，w1=0.3，w2=0.4，w3=0.2，w4=0.1，hmin=20 m，hmax=100 m，αmax=45°，βmax=45°，所有威胁区的威胁系数ρk均取1．基于文献[5-7]，BSO，QBSO，GBSO和改进IQBSO中与新个体的产生等相关参数设置为：pc=0.2，p1=0.8，p2=0.4，p4=0.5，K=5．基于文献[6]，GBSO算法的全局最优影响系数的边界值Cmax和Cmin分别设为0.8和0.2，变异系数ξ动态变化并且与搜索空间的宽度成正比，比例系数取值为0.05．基于文献[7]，在IBSO算法中，采用局部优化策略的概率pone=0.4，类别数K=2．在IQBSO算法中，随机选取两个类的概率p3=0.6，两个种群共同进化的代数Tx=150．五角星和菱形符号分别表示起始点位置和目标点位置．分别使用上述5种算法在3个不同的测试样例上进行实验．每种方法分别独立运行50次．其中，第1个测试样例对应的战场环境中威胁区较少，第2，3个测试样例对应的战场环境中威胁区分布较密集．将3个样例中全局较优航迹的代价区间分别定义为[770，1100]，[900，1200]及[860，1300]．以每种算法在50次实验中规划出的最优航迹的总代价位于全局较优代价区间的次数作为衡量算法，搜索出较优航迹的能力．5种算法在3个样例上的最优航迹的总代价的均值、标准差、最大值、最小值和较优航迹次数(W)如表1~3所示．10.13245/j.hust.240416.T001表15种算法在样例1上的最优航迹的总代价比较算法总代价W均值标准差最大值最小值IQBSO9.340×1028.610×1011.253×1037.788×10249QBSO9.959×1021.248×1021.386×1038.234×10241BSO2.259×1041.618×1046.178×1041.700×1030IBSO1.120×1032.758×1021.967×1037.809×10232GBSO1.134×1041.149×1044.413×1041.023×103210.13245/j.hust.240416.T002表25种算法在样例2上的最优航迹的总代价比较算法总代价W均值标准差最大值最小值IQBSO9.399×1022.777×1011.088×1039.265×10250QBSO1.012×1031.425×1021.418×1039.267×10242BSO3.637×1035.531×1032.282×1041.184×1031IBSO2.084×1034.025×1032.153×1049.281×10225GBSO1.856×1032.803×1032.121×1049.606×102510.13245/j.hust.240416.T003表35种算法在样例3上的最优航迹的总代价比较算法总代价W均值标准差最大值最小值IQBSO1.017×1031.463×1021.467×1038.672×10248QBSO1.161×1032.430×1021.925×1038.671×10240BSO4.993×1041.111×1054.018×1051.298×1031IBSO1.269×1033.703×1022.641×1038.761×10234GBSO1.444×1031.816×1021.785×1031.113×10315从表1~3可以看出：IQBSO算法规划出的最优航迹总代价的均值和标准差均最小，具有最高的收敛精度和稳定性．此外，IQBSO算法在3个样例上生成较优航迹的次数均最多，具有最高的全局优化能力和稳定性．BSO算法所规划出的航迹总代价的均值最高，全局搜索能力较差．IBSO和GBSO算法通过运用全局最优个体引导等方式，在一定程度上提升了BSO算法的优化能力．QBSO方法规划出的航迹总代价的均值和标准差相对较小，说明QBSO算法通过将量子机制引入到BSO算法，提升了BSO算法的全局搜索能力和稳定性．IQBSO算法通过将协同进化机制和所提出的待变异个体产生方式引入到QBSO算法，进一步提升了全局搜索能力、收敛精度和稳定性．图4所示的是5种算法在3个样例上平均最优代价值(Q)的收敛曲线，n为计算次数．可以看出：BSO，IBSO及GBSO算法在进化前期收敛速度很快，但很容易陷入局部极值；QBSO算法的全局优10.13245/j.hust.240416.F004图45种算法在3个样例上平均最优代价值的收敛曲线化能力明显优于上述3种算法，但在迭代过程中仍有可能陷入局部极值．IQBSO算法可进一步提升整个种群的多样性，从而提升了全局搜索能力．图5所示的是5种算法在3个样例中生成的最优航迹的平面投影．可以看出：IQBSO算法规划出的航迹可以很好地实现威胁规避；BSO和GBSO算法规划出的航迹很容易穿过威胁且经常违背UAV的机动约束，其总代价的均值很高．10.13245/j.hust.240416.F005图55种算法在3个样例中生成的最优航迹的平面投影图6所示的是IQBSO算法在3个样例上生成的最优航迹的剖面，图中：l为水平距离；h为高度．由图6可以看出：IQBSO算法规划出的航迹可以很好地实现地形跟随，从而提高UAV的生存概率和10.13245/j.hust.240416.F006图6IQBSO算法在3个样例上生成的最优航迹的剖面突防概率．4 结语针对QBSO算法存在的早熟收敛等问题，提出了一种改进的QBSO算法并将其应用于UAV三维航迹规划．与BSO，QBSO，IBSO和GBSO算法相比，改进QBSO算法具有更高的全局搜索能力、收敛精度和稳定性．改进QBSO算法在进化过程中仍有可能陷入局部极值．如何提高所提改进QBSO算法的优化能力并验证其优势须进一步研究．可以将所提QBSO算法与深度强化学习等人工智能算法结合，运用深度强化学习的感知能力和决策能力提高所提改进QBSO算法的优化能力．此外，改进QBSO算法在种群分类过程中采用了K-means算法，然而K-means算法的计算复杂度较高．如何降低分类算法的复杂度并提高分类效果值得进一步研究．可以将随机分类方法与基于适应度值的分类方法等进行结合，根据不同分类方法的特点在算法的不同阶段采用不同的分类方法，从而降低所提算法的复杂度．改进QBSO算法参数的值是通过参考关于传统BSO及其相关改进算法的研究成果，经过反复试验设置的．未来可以运用深度强化学习、大数据技术等对算法参数值进行优化，从而提升改进QBSO的优化能力．