智能化时代无人机产业作为迅猛发展的新兴战略产业，在货物运输、农林植保、设备监控和气象监测等方面展现出巨大的优势．随着电子商务发展、配送需求加剧，传统地面配送模式逐渐显现出弊端，尤其在“最后一公里”的城区终端配送中造成的道路拥堵、空气污染和低配送效率等问题不容忽视，亟需探索低空空域发展潜力[1]．无人机物流作为城市空中交通的重要组成部分，因其安全智能、机动灵活、绿色环保等优势将逐渐成为前景广、受众大的新型配送模式[2]．因此，研究城市复杂低空环境下物流无人机航路规划方法对于驱动无人机技术发展，促进无人机物流产业融入城市空域具有重要意义．航路规划可转化为状态空间的路径搜索问题，目前常见的有A*算法[3]、快速扩展随机树(RRT)算法[4]和智能仿生算法[5]等，国内外均已开展深入研究．文献[6]考虑城市低空空域无人机运行风险，采用改进蚁群算法规划无人机航路并进行仿真验证；文献[7]以最小化路径长度、能耗和危险度为目标，提出低空复杂环境下基于A*算法的三维无人机路径快速搜索策略；文献[8]通过建立混合整数线性规划模型，利用遗传算法规划无人机路径；文献[9]通过粒子群算法对多架无人机规划出最优防撞轨迹；文献[10]提出多约束物流无人机路径规划模型并通过B样条法对路径进行平滑处理．以上研究多集中在对航路距离代价和平滑度的优化上，未考虑无人机的航路运行风险，且规划时考虑因素较单一，未全面考虑无人机自身性能和外部运行环境的约束；A*算法和蚁群算法则对于大规模路径搜索问题收敛速度较慢，遗传算法初始状态随机性较强，可能陷入局部最优．元胞自动机作为时间、空间均离散的动力系统，通过元胞演化模拟复杂低空城市环境，根据步长有向搜索路径，提高了搜索效率．然而，传统的元胞自动机算法通常以路径长度最短为目标，并未考虑复杂低空空域限制、无人机运行风险和性能约束．本研究基于城市低空三维环境，以保障配送航路的安全性和高效性为目标，提出一种改进元胞自动机算法，以获得满足安全与距离代价双重约束的物流无人机配送航路．1 航路风险建模1.1　任务假设针对某城市区域，已知无人机物流配送点和需求点位置，现须在满足安全的前提下，规划一条距离短、风险小的航路供无人机航行．假设无人机为可垂直起降的多旋翼无人机，在飞行过程中采用分层航路，即同一条航路除了在起降点位置上无人机须要上升或下降高度到指定航路高度层，航路其余位置高度保持不变，如图1所示．10.13245/j.hust.238475.F001图1无人机分层飞行示意图1.2　运行环境构建无人机物流任务所处空间是以ABCD-A'B'C'D'为顶点，长度为L、宽度为W、高度为H的立方体．将矩形用粒度为lg的元栅格进行分割，可将任务空间划分为u×v×w个栅格，其中u=int(L/lg),v=int(W/lg),w=int(H/lg)，int(∙)为取整函数．每个栅格的中心视作备选航路点．在无人机运行空域内，存在障碍物、通信信号盲区、低空强对流气候区等多种特殊区域，威胁无人机正常运行．为了提升航路安全裕度，划设特殊区域保护区．保护区范围为特殊区域向外扩展wu/2所形成的区域(wu为无人机轴距)，无人机在执行任务过程中禁止飞入．为定量描述空域栅格特征，根据栅格地理信息特点设置栅格属性值C，如表1所示．10.13245/j.hust.238475.T001表1栅格属性值类别飞行准则C政策性禁飞区禁止飞入1通信信号盲区禁止飞入1低空强对流气候区禁止飞入1地形、障碍物实体禁止飞入1自由空域可飞入0城市复杂低空空域中，由于无人机执行任务各区域建筑分布、气象环境、地形地势等因素不同，空域栅格属性值存在差异，因此栅格危险度[7]也不同，栅格危险度rair定义为rair=1(保护区栅格);1k∑i=1kCi(可行栅格),式中k为当前栅格的邻域栅格总数．1.3　运行风险评估无人机运行风险主要体现在由于无人机机体、动力装置、飞控系统等子系统失效导致无人机坠地所造成的人员伤亡等事故．无人机运行风险表示为发生地面撞击事故的可能性与严重度的乘积，可能性体现为因无人机撞击所造成的地面人员伤亡概率，严重度体现为因撞击事故受到影响的人员数量Ncas，Ncas=PcasNeff式中：Pcas为撞击事故所致的人员伤亡概率；Neff为受撞击事故影响的人员数量．Pcas与无人机重力势能、撞击动能、下坠阻力和地面遮蔽系数有关[11]，Pcas=1/[1+μ/λλ/Ec1/(4Fs)],式中：Fs为遮蔽系数；μ为Fs=0.5时伤亡率达到50%所需的撞击动能，μ=1×106；λ为Fs趋近于0时人员伤亡的撞击动能阈值，λ=100；Ec为撞击动能；Fs∈[0,1]为地面遮蔽系数，具体取值参见文献[12]．用无人机最大运行速度的1.4倍近似替代无人机的撞击速度[12]，则撞击动能可近似表示为Ec=mvu2，式中：m为无人机最大起飞质量；vu为无人机的最大运行速度．受撞击事故影响的人员数量等价于无人机撞击事故影响面积与区域人口密度的乘积，Neff=Seffρ，式中：Seff为无人机撞击事故的影响面积；ρ为当前区域的人口密度．无人机坠地影响面积与地面人员、无人机参数和下坠倾角有关(如图2所示)，Seff=π(wu/2+rp)2+hp(wu+2rp)cotφ，式中：rp为人体半径；hp为人体高度；φ为无人机下坠倾角．10.13245/j.hust.238475.F002图2无人机撞击事故影响区域由于规划区域内地面环境的差异，因此不同栅格的Ncas数量级可能存在较大差异，须对其进行归一化处理，归一化后的无人机运行风险表示为rgro=Ncas-min(Ncas)max(Ncas)-min(Ncas)，式中min(∙)和max(∙)分别为取最小值和最大值函数．2 改进的元胞自动机算法2.1　传统元胞自动机算法原理元胞自动机是为了模拟生命系统的自复制功能所创立，由元胞、空间、状态、邻居和演化规则组成，可表示为四元组A(Cq,S,N,R)，式中：A为元胞自动机系统；Cq为q维空间；S为元胞有限的离散状态；N为邻域元胞集合；R为演化规则．根据元胞邻域结构的不同，邻域形式可分为冯∙诺依曼(von∙Neumann)型、摩尔(Moore)型和扩展摩尔型，如图3所示．10.13245/j.hust.238475.F003图3元胞邻域结构2.2　改进元胞自动机算法设计方案2.2.1　步长演化规则在航路搜索过程中，每个栅格可视作一个元胞．依据规划区域地理信息，赋予元胞初始状态值：起点元胞状态值为2，终点元胞状态值为3，保护区栅格状态值为-1，其余元胞状态值与栅格属性值相同．设航路起点和终点元胞坐标分别为(ista,jsta)和(iend,jend)，按照如下规则对元胞进行演化：从规划区域起点元胞开始，向终点元胞进行演化．设置规划区域起点元胞为当前元胞，当前元胞的状态值遵循演化公式Sij=Cij    (Cij≠0);min(CMoore(i,j)≥2)+1    (Cij=0),式中：Sij为当前元胞i,j的状态值；CMoore(i,j)为当前元胞摩尔邻域内可行元胞的状态值集合．重复以上演化过程，当终点元胞邻居中所有可行元胞的状态值全部发生变化时演化终止．图4所示为元胞演化示意图，演化终止的判别公式为CMoore(iend,jend)=0=∅．10.13245/j.hust.238475.F004图4元胞状态值演化2.2.2　航路规划约束a. 转弯角约束．由于无人机最大转弯半径的限制，须对航路转弯角进行约束．航路转弯角φi须满足0≤φi≤φmax，φmax为无人机最大转弯角．b. 转弯缓冲距离约束．无人机在转弯之前需要一定的缓冲距离，才能确保安全过弯．假设航路上相邻的两个转弯点坐标为(xi,yi,zi),(xi+1，yi+1,zi+1)，默认在转弯点之间沿直线飞行，规划航路时须满足[(xi-xi+1)2+(yi-yi+1)2+(zi-zi+1)2]1/2dmin，式中dmin为无人机转弯最小缓冲距离．c. 步长连续性约束．可行元胞状态值减去2即可得到元胞步长．当航路回溯搜索时，两个连续航路点间的元胞步长Tij和Ti+1,j+1须连续而不能跳步，表示为Tij-Ti+1,j+1=1．2.2.3　初始航路回溯搜索获得全局可行元胞步长后由终点向起点回溯搜索初始航路．安全航路以物流无人机安全运行为目的，在满足环境约束和任务要求的前提下提高飞行效率，缩短航路长度．航路风险体现在航路危险度与无人机运行风险共同作用的结果，因此将航路风险值看作航路危险度与无人机运行风险值的和．寻路时引入代价函数P作为选择元胞的依据，代价函数可表示为P=α1D+α2R，式中：D=d1+d2为航路距离代价，d1为起点元胞到备选元胞间的欧氏距离，d2为终点元胞与备选元胞间的欧氏距离；R=rair+rgro为航路风险代价； α1,α2分别为航路长度和风险的启发因子．在回溯搜索选择航路点时应尽可能选择代价函数值最小的元胞，表示为Pi+1,j+1=min(PMoore(i,j))，式中：PMoore(i,j)为当前元胞摩尔邻域内满足航路规划约束的可行元胞代价函数值集合；Pi+1,j+1为下一步元胞的代价函数值．当Tij=1时航路回溯搜索结束，得到初始航路点集合W1,W2,…．2.2.4　航路优化栅格粒度影响了搜索效率和航路精度．栅格粒度越大，航路搜索效率越高，但航路精度会有所下降，反之亦然．为在提升搜索效率的前提下不以牺牲航路精度为代价，对算法搜索的初始航路进行优化，具体思路为：对于初始航路中的一系列航路点Wi,Wi+1,…，若航路点Wi,Wi+2之间的航路不存在保护区，则移除航路点Wi+1，将航路点Wi,Wi+2直接相连；判断新航段Wi-Wi+2与航段Wi+2-Wi+3之间是否满足航路转弯角约束，若满足则保留优化航段，否则恢复原始航段．重复以上步骤，判断航路点Wi与Wi+3,Wi+4,…之间的航路是否存在保护区，若存在则终止本次迭代，依次对Wi后的航路点执行以上算法直至航路终点；否则不断迭代移除中间航路点直至到达航路终点．图5所示为航路优化算法示意图，若初始航路点依次为W1,W2,W3,W4,W5,W6，则经过优化后的航路点为W1,W3,W4,W6．10.13245/j.hust.238475.F005图5航路优化过程2.2.5　算法实现改进元胞自动机算法分三部分实现：首先，基于城市低空环境特征更新元胞状态值，根据规划区域保护区数据和起始点位置计算全局可行元胞的步长；然后，根据全局可行元胞步长和代价函数值搜索最佳航路点，生成初始航路；最后，采用航路优化算法对初始航路进行优化从而得到最终航路．算法具体流程如图6所示．10.13245/j.hust.238475.F006图6算法流程图3 仿真与分析3.1　仿真实验环境为验证本研究提出的无人机航路规划算法的合理性和有效性，选择南京航空航天大学将军路校区作为无人机航路规划背景．通过实地测绘获取该区域1 025 m×1 060 m×70 m的地理信息数据，采用栅格法对规划环境建模，人口密度ρ取0.025人/m2，无人机轴距wu取1 000 mm，无人机最大运行速度vu取18 m/s，无人机质量m取15 kg，人体半径rp和高度hp分别取0.3 m和1.7 m，无人机最大转弯角φmax取90∘，无人机转弯最小缓冲距离dmin取5 m，航路长度和风险的启发因子α1和α2均取0.5．3.2　仿真结果分析使用Matlab进行仿真实验来验证本文算法的有效性．在区域内选取3对配送路线，起点和终点分别记为O1,O2,O3和D1,D2,D3．按照改进算法规划的无人机航路如图7所示．10.13245/j.hust.238475.F007图7航路规划结果为进一步验证改进元胞自动机算法的优越性，分别采用传统元胞自动机算法、改进元胞自动机算法和A*算法对航路进行规划，规划的航路结果如图8所示，航路参数见表2．10.13245/j.hust.238475.F008图8传统算法与改进算法对比10.13245/j.hust.238475.T002表2各规划航路主要参数对比参数航路传统算法改进算法A*算法航路长度/kmO1-D11.306 71.122 11.170 0O2-D21.277 21.181 11.220 5O3-D31.015 00.846 40.870 6航路风险O1-D11.982 51.898 82.135 0O2-D21.910 01.831 12.017 5O3-D31.747 51.163 81.713 8航路代价O1-D11.644 61.510 41.652 5O2-D21.593 61.506 11.619 0O3-D31.381 31.005 11.292 2由表2可知：三种算法均可完成航路规划任务，相比传统算法；改进算法规划的3条航路长度平均减少12.76%，航路风险平均减少13.92%，航路代价平均减少13.63%；相比A*算法，改进算法规划的3条航路长度平均减少3.37%，航路风险平均减少17.46%，航路代价平均减少12.6%．因此改进算法在减少航路距离的同时，提高了运行安全和效率．由图8可知：传统算法规划的航路较为曲折，且转弯点较多，而改进算法中的优化算法改善了航路结构，减少了不必要的转弯点．将考虑风险与未考虑风险的规划航路进行对比，即航路长度启发因子α1取值分别为0.5和1.0，结果见表3．10.13245/j.hust.238475.T003表3航路规划算法性能对比类别航路航路长度/km航路风险考虑风险O1-D11.122 11.898 8O2-D21.181 11.831 1O3-D30.846 41.163 8未考虑风险O1-D11.032 82.142 5O2-D21.147 21.857 0O3-D30.787 51.230 0以O1飞行至D1的航路为例，考虑航路风险的规划航路长度为1.122 1 km，比未考虑航路风险的航路长度增大0.089 4 km；航路风险从2.142 5降低至1.898 8，降低比例约11.38%．考虑航路风险的规划航路相比未考虑航路风险的规划航路在长度上略有增长，但航路风险有明显下降，有利于物流无人机在复杂城市低空环境下的安全运行．4 结语为保障物流无人机在城市低空复杂空域的安全、有序运行，提出一种基于代价函数的改进元胞自动机航路规划算法．该算法以航路长度和风险的线性组合作为代价函数，相比传统算法和A*算法在规划航路的长度和风险上优势显著．仿真实验表明：航路长度平均减少12.76%，航路风险平均减少13.92%；考虑航路风险的改进算法相比未考虑航路风险的算法，航路长度略有增大，但航路风险有较大下降，使得航路在距离代价增长较小的条件下，较大程度提升航路的运行安全．