复杂低空物流无人机路径规划

张启钱，许卫卫，张洪海，邹依原，陈雨童

（南京航空航天大学 民航学院，南京211106）

无人机作为科技创新的重要产业处于井喷式发展时期。由于技术日臻成熟，各型别无人机应用领域不断扩大，如军事侦查、农林生产、物流运输［1］。虽然其在物流领域的应用处于起步阶段，但各国进行了前沿性研究且成功案例较多。国内外学者对军用无人机研究深入［2-3］，涉及民用无人机则较少，且主要对巡检、药物喷洒等路径进行规划［4-5］。国内学者研究物流无人机管控居多［6］，国外学者则关注其路径规划，解决如何将货物送至多个分快递站问题［7-8］，但却很少研究分站点间具体运输的路径。实际应用中，物流无人机自身及环境复杂，其路径是货物安全、经济、快捷送至目的地的重要保障，可降低运营成本、提升流转效率，故研究物流无人机路径规划技术有重要现实意义。

现有对复杂低空物流无人机路径规划的研究较少。Byung等［7］融合无人机物流系统飞行时间有限、货物对性能影响大等特点，提出基于混合整数线性规划的无人机运输路径模型，设计启发式算法解算路径；Boualem等［8］考虑无人机在人道救援物流方面的应用，提出在载荷和能耗约束下无人机运输物品总成本最小的优化模型；Scott等［9］验证了无人机运输医疗物资可行性，提出2个面向医疗服务的无人机物流配送模型；周浪 提出“配送车+无人机”配送及路径优化问题，构建了不同运输路径优化模型，采用遗传算法等求解模型；翁丹宁［11］、吴永鑫［12］探讨无人机物流配送的短板及未广泛用于市场的原因，分析物流无人机的优势并研究了影响其配送的因素；Jiang等［13］考虑无人机性能建立带时间窗的无人机物流任务分配模型，设计改进粒子群算法进行求解；Dorling等［14］推导并证明多旋翼无人机能耗与载荷呈线性变化，建立无人机路径规划的混合整数线性模型；Torabbeigi等［15］构建并验证无人机电量消耗率随载重的变化函数，采用最小集覆盖法对物流无人机路径规划进行建模，提出一种变预处理算法以提高求解模型的速度；Goss等［16］建立考虑无人机物理特性和执飞任务特点的无人机性能模型，采用四旋翼无人机测试并验证模型的正确性；Abeywickrama等［17］研究无人机机动行为的电量消耗，考虑不同场景及条件对无人机能耗的影响，提出涵盖无人机飞行全过程的能耗模型。

纵观已有成果，部分将物流无人机路径规划问题简化为道路车辆路径规划问题，未体现无人机空中运动特点及性能约束［12-15］；部分在规划物流无人机路径时未考虑实际飞行环境且未将任务要求、能耗、货物质量等影响要素组合考虑，规划结果不理想［7-8，10］。现有研究面向飞行行为、宏观路径规划，考虑的路径影响要素单一有限；而低空环境下物流无人机路径规划空间复杂，路径搜索需结合外部环境与自身限制，且对搜索效率要求高。同时物流无人机路径规划与执行其他任务的路径规划不同，其航空物流特点导致规划考虑的因素更全面、具体，如何将这些因素综合考虑以进行路径规划亟待研究。

本文从物流无人机飞行安全和运输效率角度出发，以最小化路径飞行时间、能耗及危险度作为目标函数，构建在满足无人机运输货物过程中自主安全避障的前提下，减小运行成本的物流无人机路径规划模型，设计基于A*算法的路径快速搜索策略，以获得最佳的货物配送路径。

1 物流无人机路径规划问题建模

1.1 问题描述与相关假设

某区域有多个分快递站且都配备物流无人机，利用无人机完成快递站间货物运输任务。假设物流无人机均为可垂直起降的充电旋翼式无人机，且对货物载重、能耗等有限制。为确保货物安全快捷运输至分快递站，要求对其路径进行科学规划以降低运输风险及成本。已知物流无人机有2种运输模式：①每次为单个分快递站（即目的派送点）提供服务后返回配送中心（即起始派送点）；②为多个分快递站提供服务后再返回配送中心。本文采取第1种模式研究单架物流无人机从起始派送点到目的派送点运输代价最小的安全避障路径规划技术。

1.2 飞行环境建模

常采用函数模拟法和地形高程数据对地形地物建模。前者通过解析公式对实际地形进行计算机模拟，如文献［18］中的函数，但该方法太过理想化，利用函数产生的地形与真实地形存在差距，这对于以安全为主的物流运输场景不适用；后者通过对已有地理信息数据进行插值处理获得逼近真实环境的地形。除避开地形地物外，各地政府部门对无人机的活动范围也进行限定，规定了一系列无人机禁飞区、限飞区和危险区等［19］。

1.3 栅格法规划环境表征

栅格法将环境划分为单元格，依据其中是否存在障碍物分类：若存在地形地物等障碍，则为障碍栅格并赋值1；否则，为自由栅格并赋值0。采用栅格法扩展路径点需考虑栅格粒度大小设置，即栅格边长。栅格粒度过大导致规划的路径偏离理论最优解，过小则导致规划时计算量大。由于栅格中心点是物流无人机飞行路径点且路径需满足无人机性能约束，因此栅格粒度大小需和性能限制匹配。

设运输任务规划环境是长、宽、高分别为x、y、z的长方体区域，记为OABC-O′A′B′C′，以O为原点建立三维直角坐标系。其中，OA长度为x，AA′长度为y，OC长度为z。考虑环境信息及物流无人机性能确定栅格粒度大小l。按以下2步确定栅格坐标。

步骤1 对OABC-O′A′B′C′沿OO′进行m等分，过等分点作平行于OABC面的平面，得m-1个平面γj（j＝1，2，…，m）。其中，m ＝int（y／l），int（）为取整函数。

步骤2 对上述任一平面γj，沿OA进行u等分，沿OC进行h等分。其中u＝int（x／l），h＝int（z／l）。此时γj便离散成u×h个平面栅格，OABC-O′A′B′C′被离散成m×u×h个立体栅格，则物流无人机待扩展路径点可标记为p（i，j，k），i＝1，2，…，m，j＝1，2，…，u，k＝1，2，…，h。

图1为规划环境栅格化示意图。

图1 规划环境栅格化Fig.1 Grid of planning environment

1.4 多限制条件物流无人机路径规划模型

物流无人机货物运输需考虑多方面影响要素，本文综合考虑物流无人机性能约束与任务要求，从运输的安全性、经济性和快捷性等角度建立多限制条件物流无人机路径规划模型。

1.4.1 性能约束

1）最远航程、最小路径段长度

物流无人机路径中相邻两飞行路径点间距离不能小于最小路径段长度，约束为

式中：tfly、D分别为飞行时间、路径危险度；a1、a2、a3分别为飞行时间、能耗、危险度的系数；T为完成任务可用总时间；t1为起始派送点起飞时刻；t2为到达目的派送点最晚时刻。

式（8）为目标函数；式（9）中前7个式子为物流无人机性能约束；最后一个式子为物流无人机货物运输任务约束，表示物流无人机应在规定时间内完成任务。

2 算法设计

2.1 A*算法原理及流程

A*算法通过估计代价函数f（n）选择待扩展路径点n，通常考虑的代价是航程，即要求获得飞行总航程最短的路径，表达式为

式中：g（n）为起始派送点S至待扩展路径点n的实际飞行航程，称为实际代价函数；h（n）为待扩展路径点n至目的派送点G的预估飞行航程，称为启发式函数。

A*算法的原理及实施流程参见文献［20］，本文不再赘述。

2.2 算法设计方案

直接套用A*算法求解上述模型时存在以下问题：①函数g（n）及h（n）仅计算航程，却没考虑能耗、飞行时间等因素，不符合实际应用；②忽略自由栅格位置的差异而同等看待，低空环境及障碍物对路径点选择的影响未体现；③未考虑无人机物理性能、机动行为对路径代价的影响；④函数f（n）未体现无人机货物载重；⑤待扩展路径点时前期和后期的效率、精度不匹配；⑥规划的路径存在冗余路径点且不平滑。

为适用本文模型，采用以下方案设计求解算法：①物流无人机飞行航程相比于车辆运输距离已大大缩短，因此航程不是计算的重点。物流无人机飞行由电池供电且目的派送点接收货物有时间限制，所以对飞行时间及能耗等代价进行考虑。此外，爬升时高度差导致的能耗与平飞相同距离能耗不同［16］，故不同操作的成本代价也应分开计算。②低空安全飞行是物流无人机运输货物的首要目标，为有效躲避地形地物、恶劣气象等障碍，定义栅格危险度因子概念，用新方式对数字栅格环境存储，并在g（n）中新增栅格危险度因子以确保扩展更安全、稳定的路径点。③考虑所载货物对物流无人机性能、路径有影响，定义货物质量惩罚系数，对飞行时间、能耗等代价进行惩罚。④h（n）采用执行任务时间和电量表示，将其表达为从当前路径点预计到达目的派送点的时间和能耗与剩余可飞时间和电量的差值，该值大小反映完成运输任务成本的高低。⑤A*算法中，当待扩展路径点靠近目的派送点时h（n）变小，导致其在f（n）中所占比例减小，启发效率降低，此时可对h（n）适当增加比重，但这导致前期搜索范围小而无法搜索最优路径。因此，设计动态权重系数对g（n）和h（n）加权以权衡算法搜索效率和精度。⑥为筛除原路径中冗余路径点，采用双向交叉判断法对路径点对连线构成的路径段和障碍物进行相交检测以获得精简路径，并利用样条曲线插值对精简路径处理得到物流无人机连续平滑路径。

2.3 具体设计

2.3.1 新增栅格危险度因子

A*算法利用0-1矩阵对栅格环境存储，这将自由栅格同等看待，在搜索路径点时这些区域被选中的概率相等，导致搜索的路径安全性差。实际情况中，虽然某些栅格不对物流无人机造成危害，但由于栅格所处位置不同造成无人机在不同栅格中飞行所受危险度有差异。本文对原标记为1的栅格不做改变，而对标记为0的栅格将其值重新标为危险度因子。定义栅格危险度因子来衡量路径点周围环境对该点的危险程度，公式如下：

式中：d为栅格危险度因子；Nobstacle为该栅格周围障碍栅格数；Nsurround为周围总栅格数。

本文称对自由栅格差异化而非同一化处理的方法为非二值存储法，此法能保留原始威胁信息。图2为栅格危险度因子计算示意图，黑、白栅格分别代表障碍栅格和自由栅格。式（12）为A*算法0-1矩阵存储结果，式（13）为非二值存储法存储结果。

图2 栅格危险度因子计算Fig.2 Calculation of grid risk

2.3.2 添加货物质量惩罚系数

物流无人机不同机动行为对能耗是不同的［16］，所以函数g（n）将水平和升降运动能耗分开计算且不考虑起降能耗。假设物流无人机水平运动能耗为距离的λ倍，垂直运动能耗为高度差的r倍，g（n）表达式为

2.3.5 路径优化及平滑

路径转弯多表明物流无人机姿态改变次数多，利用双向交叉判断法筛除冗余路径点可减少路径点数，其思想是：当原路径中不相邻两路径点间连线构成的路径段不与障碍物交叉时，可剔除该两点间的冗余路径点。假设采用本文算法搜索的原路径为Path，按如下步骤对该路径进行优化，操作示例如图3所示。

步骤1 创建新列表DELETE。

步骤2 对Path中路径点按其高度划分为多段子路径Path＝｛subpath1，subpath2，…｝，任一子路径可表示为subpathi＝｛p1，p2，…，pe｝，e为该段路径点总个数。

步骤3 对subpath中满足e＞2条件的任一子路径进行操作：①共需进行q轮判断，其中q＝e-2，并令i＝1；②若i≤q，则跳至步骤4判断点pj与pj+e-i是否位于DELETE中（j＝1，2，…，i），若i＞q，跳至步骤6。

步骤4 对上述每一路径点判断其是否位于DELETE列表中。若路径点pj与pj+e-i均不在DELETE中，则跳至步骤5。

步骤5 若路径点对间连线与障碍物不交叉，则该两点可直接连接，并将原段中两点间的其余点移入DELETE，第i轮中路径点对全判断完后跳至步骤3中的②，令i＝i+1。

步骤6 对其余子路径循环执行步骤3～步骤5，直至所有路径点全判断完毕。

步骤7 将位于DELETE的路径点在原Path中筛除，剩余路径点序列构成的路径即为精简路径。

图3 双向交叉判断法Fig.3 Bidirectional cross judgmentmethod

利用双向交叉判断法得到的精简路径为一条折线飞行路径，但在复杂低空中物流无人机从起始派送点至目的派送点的路径应是连续的，故对路径再进行平滑处理。采用B样条曲线对精简路径插值拟合获得可飞的连续平滑路径，具体流程见文献［21］。

2.4 算法实现

本文算法规划物流无人机路径的流程如图4所示，步骤如下：

步骤1 获得物流无人机机动性能约束，起始派送点和目的派送点坐标，货物质量惩罚系数τ（Q）。

图4 本文算法流程Fig.4 Flowchart of proposed algorithm

步骤2 依据运输场景地图确定栅格粒度大小，采用非二值存储法生成栅格危险度因子矩阵。

步骤3 建立OPEN、CLOSE列表，将起始派送点所在栅格加入OPEN列表，将障碍栅格加入CLOSE列表。

步骤4 循环执行：

4.1 搜索前一路径点周围栅格危险度因子不为1的路径点并放入OPEN列表，计算其gQ（n）、hQ（n）、W（n）、W′（n）和fQ（n）值，选择其中fQ（n）值最小的点为当前正扩展路径点并记作A。

4.2 将A点移入CLOSE列表。

4.3 判断A点周围的路径点，若其危险度为1则跳过该点，否则记该点为B并求gQ（B）、hQ（B）、W（B）、W′（B）和fQ（B）值，进 行 以 下判断：

4.3.1 若B点已在OPEN列表，则确定B点的最优父节点。记对原先B点求得的估计成本值为fQ（B′），若fQ（B）＜fQ（B′），B点的父节点为A点；否则B点的父节点不变。

4.3.2 若B点已在CLOSE列表，则跳过B点。

4.3.3 若B点既不在OPEN列表也不在CLOSE列表，则将B 点移入OPEN 列表并在OPEN列表中对估计成本值的大小按升序排序，选出fQ（n）值最小的点为下一路径点。

步骤5 当路径点到达目的派送点所在栅格或未扩展到目的派送点但OPEN列表为空时搜索结束，终止步骤4。

步骤6 从目的派送点通过各路径点的父节点反向移动回到起始派送点，便得物流无人机路径，并采用双向交叉判断法及样条曲线插值对该路径进行优化平滑。

3 仿真验证及分析

3.1 仿真环境及参数设置

为验证本文算法规划物流无人机路径的有效性，利用Google获取某区域1 300 m×1 300 m×310m的地形高程数据，并将该区域作为此次物流无人机路径规划环境，部分数据见表1，其中，第2行第2列表示（0，0）点处地形高程值（此处为34.52m），相邻单元格间的距离为5 m，即第i行第j列单元格中数据代表坐标（5（i-1），5（j-1））处的高程值。采用MATLAB将此环境以图形化展示，如图5所示。由于目前直接面向顾客的无人机末端配送技术不成熟及相应基础设施建设不完备，本文规划的物流无人机路径仅用于两货物存储仓库间。为尽量真实模拟物流无人机最优飞行路径产生的过程，仿真实验具体参数按表2进行设置。

表1 环境数据Tab le 1 Environm en tal data

图5 物流无人机路径规划环境Fig.5 Environment for path planning of logistics UAV

表2 仿真参数Table 2 Sim ulation param eters

3.2 仿真结果分析

采用本文算法及表2中参数设置得到物流无人机路径规划结果，如图6中绿色曲线所示。可知，该物流无人机安全平滑地从起始派送点飞行至目的派送点，且有效地对低空环境中的障碍物进行避让。为对比本文算法与A*算法规划性能的优劣，在上述各参数设置、路径规划环境保持不变的前提下，对该物流无人机路径利用A*算法进行规划，结果如图6中红色曲线所示，并从飞行路径的路径点数、航程、能耗等方面进行分析，数据如表3所示。可以看出，本文算法的各指标数据均优于A*算法，且路径点数变化率最大，降低了42.9%，表明双向交叉判断法可有效减小路径点数，降低物流无人机姿态改变次数，保证了货物运输安全。

为进一步分析本文算法的规划性能，再采用人工势场法、未优化及平滑的改进算法对物流无人机的路径进行规划，各算法规划结果如图7所示（为便于观察，图中仅展示各路径），并将对比结果记录于表4。由图7知，在图7（a）中，各算法均可规划出路径，但A*算法及人工势场法规划的路径高度变换次数多，且人工势场法陷于局部最优；在图7（b）中，未优化及平滑的改进算法规划的路径存在冗余路径点且不平滑。再结合表3、表4分析，人工势场法在复杂环境下陷入局部最优解，故规划时间、路径点数、路径航程分别为本文算法的6.5、1.6、1.2倍；A*算法未考虑物流无人机货物运输要求导致其路径长、能耗多，飞行时间约为本文算法的1.1倍；未优化及平滑的改进算法在规划时间上比本文算法短，原因在于双向交叉判断法、样条曲线插值都具有一定时间复杂度，但其路径点数、路径航程、危险度均比本文算法高，故规划时间变长属可接受范围。

图6 物流无人机路径规划结果Fig.6 Path planning results of logistics UAV

表3 A*算法和本文算法路径规划结果对比Tab le 3 Com parison of path p lanning resu lts betw een A* algorithm and proposed algorithm

通过以上分析得，本文算法路径规划结果在路径航程、能耗、危险度等方面都优于A*算法及人工势场法；物流无人机路径规划模型可优化无人机的货物运输时间和路径搜索结果，实现时间上飞行时间减少、空间上安全平滑避障，进而及时低风险地完成任务，这也表明引入栅格危险度因子、动态赋权估计成本函数等改进的合理性。

图7 不同算法路径规划结果Fig.7 Path planning results of different algorithms

表4 不同算法路径规划结果对比Tab le 4 Com parison of path p lanning resu lts am ong different algorithm s

3.3 算法参数分析

在求解物流无人机路径规划模型时，栅格粒度大小l和代价权重值组合｛α1，α2，α3｝会对模型解算结果产生影响。本文采用对照实验法分析参数l和｛α1，α2，α3｝的取值对路径规划结果的影响。

保持代价权重值组合如表2中设置不变，栅格粒度大小l分别取5、10、15、20m，在其他参数设置及规划环境相同的条件下进行多组对照实验，对各组实验规划出的路径数据记录于表5中，并画出各数据的变化趋势图，如图8所示。

表5 栅格粒度大小对路径规划结果的影响Table 5 Influence of grid length on path p lanning resu lts

图8 栅格粒度大小的影响Fig.8 Influence of grid length

从表5得出，l＝5m时，路径点数为129，栅格危险度因子为11.69；l＝15 m 时，路径点数为128，与l取5m时接近，但其栅格危险度因子却比前者高65%；l分别为10m、20m时，路径点数、栅格危险度因子均明显高于l取5m情况。由图8知，图8（a）中，路径点数、路径航程随栅格粒度大小呈先增大后减小变化，但仍有上升趋势；图8（b）中，能耗随栅格粒度大小有上升变化趋势，危险度则相反。经分析，栅格粒度太小变化时，l越小，环境划设精细，待扩展路径点多，需要更多存储空间；l越大，环境划设粗糙，效率提高，但难保证航程等成本小。综上所述，l取5m时可规划出结果最佳的飞行路径，因此选择性能较好的5 m为最优的栅格粒度大小设置。

同理，保持栅格粒度大小l＝5m设置不变，由于物流无人机路径对安全要求较高，故固定α3为0.5不变，分析α1、α2取值对结果的影响，代价权重值组合取表6中设置，在其他参数设置及规划环境相同的条件下再进行多组对照实验，对各组实验规划出的路径数据记录于表6中，并画出其变化趋势图，如图9所示。

由图9知，各路径数据随代价权重值的取值不同有明显变化趋势。图9（a）中，路径点数虽有起伏但呈下降趋势，路径航程变化幅度大且在实验5时取最小值；图9（b）中，能耗呈先减小后增大趋势，且在实验5后增长率达7.3%，危险度呈下降变化，且在实验5中取相对较小值。经分析，飞行时间权重α1越大，路径点数相应变少，但能耗等有所增加；能耗权重α2越大，能耗相应变少，但路径点数、危险度有增加趋势。综上，实验5中参数设置可平衡各路径代价，故选择α1＝0.4，α2＝0.1，α3＝0.5为最优的代价权重值组合。

表6 代价权重值对规划结果的影响Tab le 6 In fluence of cost weight on p lanning resu lts

图9 代价权重值的影响Fig.9 Influence of cost weight

4 结 论

1）引入物流无人机物理性能及任务约束，设计了考虑飞行时间、能耗及危险度的目标函数，建立多限制条件物流无人机路径规划模型；为快速解算该模型，设计A*算法并采用双向交叉判断法、样条曲线插值对原路径优化平滑。

2）以某区域地形高程数据进行仿真验证，结果表明，本文算法求解物流无人机路径规划模型时能在无人机安全避障的前提下，缩短规划时间，减少能耗与路径点数，提供可飞连续平滑的路径。

3）本文模型及算法可应用于物流无人机货物运输路径的规划，且规划结果受栅格粒度大小和代价权重值影响较大。在本文设置的运输路径规划环境下，当栅格粒度大小取5 m，代价权重值取0.4、0.1、0.5时，路径规划结果最佳。

对低空物流无人机路径规划技术的研究目前较少，下一步的研究将从突发情况下物流无人机实时动态路径规划及多物流无人机协同货物运输路径规划进行。