复杂边界田块旋翼无人机自主作业路径规划

黄小毛 张 垒 TANG Lie 唐 灿 李小霞 贺小伟

(1.华中农业大学工学院，武汉 430070；2.爱荷华州立大学农业与生物系统工程系，埃姆斯 IA 50011；3.华中农业大学信息学院，武汉 430070；4.塔里木大学现代农业工程重点实验室，阿拉尔 843300)

0 引言

与地面机器相比，无人机空中作业具有通过性强、效率高、精度高等优点，特别适合于植保、遥感、巡查和播撒等作业环节。作业路径是无人机自主作业时机器轨迹控制的跟踪目标，在很大程度上决定机器作业过程的作业质量、效率和总消耗。因此无人机作业路径的规划和优化工作具有重要意义。

在目前农机领域作业路径的相关研究中，针对地面作业机器的居多[1-18]，针对农用无人机作业路径的研究不多。MOON等[12]针对喷粉无人直升机，提出基于田块分区处理操作的路径规划方法。BARRIENTOS等[13]将待作业区域进行网格化处理，提出一种基于任务协商机制的多架无人机协作图像采集作业路径规划算法。VALENTE等[14]在此基础上提出一种基于“和谐搜索”(Harmony search,HS)的启发式搜索算法，求解遍历网格的最优次序。国内相关的研究目前还处于起步阶段，徐博等[15-16]对五边形边界单田块的无人机作业路径提出了具体的算法，并对多田块路径调度优化进行了研究。徐东甫等[17]、刘浩蓬等[18]分别将简单四边形边界的规划路径应用到无人机赤眼蜂投放和喷雾试验中。王宇等[19]提出一种基于栅格模型的路径规划算法。彭孝东等[20]针对固定翼无人机作业过程中转弯方式进行了详细的讨论。以上研究的共同特点是作业田块的边界较为简单，均可简化为凸多边形，而非实际应用中遇到的复杂边界田块。

本文针对旋翼农用无人机作业时可能遇到的各种复杂田块边界形状的问题，提出一种对田块边界形状具有普适意义的旋翼无人机作业路径规划算法，以获得凸多边形、凹多边形、带孔洞多边形或多个多边形的复杂边界田块下的飞行作业轨迹。

1 作业路径规划的基本原理

1.1 旋翼无人机特点及路径规划的基本要求

与大田地面作业机器及固定翼飞机相比，旋翼农用无人机具有不与地面接触、可悬停、可任意方向移动等特点，因此其作业路径规划应满足以下基本要求：有效工作路径(即作业航线)做到对待作业田块所有区域的全覆盖，重复作业及遗漏作业面积尽可能小；航线间跳转转移等非工作路径总长度尽可能小；无需考虑转弯及碾压问题；算法的时间消耗不能太长；具体作业类型及应用场合的其他要求(如避障)。

1.2 数据定义及算法流程

图1 田块边界及顶点存储顺序示意图

本文只针对常见的二维田块，用一组多边形分别代表田块的外边界和内边界，且外边界以逆时针方向、内边界以顺时针方向依次存储各顶点(如图1所示)。通过分析可知采用多边形表征田块区域时具有以下特点：每个田块有且只有一个外边界，而可能没有也可能有多个相离的内边界，且内边界一定包含在外边界内部；每个边界都是一个简单多边形(除相邻边外，其他各边之间互不相交)，可以是凸多边形也可以是凹多边形；每2个田块之间的边界均为相离关系(既不相交也不相互包含)。这些特点既是后续算法中导入的田块边界的合理性规定，也是计算机能够对多边形边界进行分组进而自动识别田块的依据[21]。

以多边形边界和作业基本参数为输入条件，以获得旋翼无人机飞行作业路径为目标，本文设计的算法基本流程如图2所示。

图2 算法流程图

1.3 边界合规性检查及分组

田块边界一般通过人工手持GPS直接打点或在GIS系统中鼠标打点获得，打点、存储的顺序可能具有随意性，不一定完全符合1.2节中的规定。因此经纬度球面坐标转换成直角坐标后，首先需要对表征边界的多边形进行合规性检查，对不合要求的地方进行调整，以保证后续处理顺利进行。然后，通过分组算法[21]确定多边形之间的关系，即从数据结构上建立田块的概念。如从图3a中输入的B0、B1、B23个边界多边形，通过分析它们之间的包含关系，分组后分别建立F0(图1)和F1(图3b)2个独立的田块区域。

图3 边界合规性检查及分组

1.4 填充线求解

多边形填充是计算机图形学中的经典问题[22]，在计算机动画、数控型腔铣削、3D打印和机器人等领域均有应用[21]。利用扫描线及多边形边线的连贯性原理，采用有序边表法进行快速求交解算，可避免传统做法中对每一组线段进行逐一相交性测试，从而大大提高了求解效率和速度。

对于给定的任意形状的单一多边形，用一组水平或垂直的扫描线进行扫描，对于每一扫描线而言均可求出与多边形边的交点，且交点成对出现，位于多边形内部的交点线段即为填充线。对于孔洞多边形而言，分别求取外边界多边形B0、内边界多边形B1～Bn的填充线，得到填充线段集合{L0}～{Ln}，再通过减法布尔运算，即可得到孔洞多边形区域的填充线集合，即航线集合{Pi}。对于孔洞多边形田块Fi而言，航向角θ条件下的作业航线Pi(θ)计算式为

(1)

以田块F0为例，外边界多边形B0与内边界多边形B1的填充线分别如图4b浅蓝色线条和图4c红色虚线条所示，二者作减法布尔运算后的结果即为田块F0在航向角为0条件下的作业航线，如图4a浅蓝色线条所示。

图4 孔洞多边形填充线的求解过程

1.5 活性边表法

单个多边形的填充线求解时，扫描线与多边形边的交点具有一定规律性，因此求解时不需要将每一条扫描线与每一条边进行求交测试并逐一解算交点，而是运用活性边表法快速求解。实施过程中用边表(Edge table，ET)来确定哪些边是下一条扫描线求交计算时应及时加入运算的边，将与当前扫描线相交的边称为活化边(Active edge list，AEL)，其个数可能不唯一且可动态变化。由活化边组成的表称为活性边表(Active edge table，AET)，基于活性边表的多边形填充线求解算法称之为活性边表法或有序边表法。

作为存储多边形各边的一种数据结构，边表用来表示对于某条扫描线而言第1次出现的边的信息。以水平扫描线(即航向角为0)为例，每个节点的信息为(Ymax，dX，XYmin，next)，其中Ymax为对应边的最大Y值，dX为沿该边从当前扫描线到下一条扫描线之间的X方向增量(当扫描线间距为1时为该边斜率的倒数)，XYmin为该边的下端点的X坐标，next为指向下一条边的指针。建立边表时，先按下端点的纵坐标(Y值)对所有边作分类和排序，再将同一组中的边按下端点X坐标递增的顺序进行排序，X坐标还相同的按dX递增的顺序进行排序。图4b中多边形的边表如表1第2列所示(表中“→”表示指针next)。

表1 活性边表法实施过程示例

对应的，活性边表AET中每个节点的信息也由4部分组成，即(Ymax，dX，XYmin，next)，其中Ymax、dX、next指针的含义同边表。当作业幅宽w不等于1时，dX=ΔX/w，式中ΔX为对应边在两扫描线间的X方向增量。与边表ET节点信息的主要差别在于第3个信息，XYmin为边与当前扫描线交点的X坐标而不是当前边的下端点X坐标。活性边表的生成是基于边表的，而由活性边表生成成对的交点即为填充线段的端点。图4b中多边形的活性边表及交点如表1第3、4列所示。

活性边表法求解多边形填充线段的流程如图5所示。对于每一条扫描线而言，对应的处理步骤如下[23]：

(1)对于扫描线Y=y，若对应的ET中非空，则将其所有的边从ET中取出并且插入到边的活性边表AET中，并对AET中各边按XYmin及dX或ΔX递增排序。

(2)若相对于当前扫描线的AET非空，则将AET中的边两两依次配对，即第1、2边为一对，第3、4边为一对，依此类推，每一对边与当前扫描线的交点所构成的区段位于多边形内，即为对应的填充线。

(3)将当前扫描线的纵坐标Y累加w，即Y=Y+w。

(4)将AET中满足Y=Ymax的边删去。

(5)将AET中剩下的每一条边的X域累加dXw，即X=X+dXw。

(6)重复执行步骤(1)，直至AET中边全部删除。

图5 活性边表法求解多边形填充线段算法流程图

1.6 航线排序优化

上文中每一填充线即为飞行器作业时的具体航线，由于填充线段的长度、位置各不相同，因此航线间存储顺序将决定遍历时航线跳转距离和遍历过程的总轨迹长度。若将每一条航线视为一个城市，这将是一个典型的旅行商问题(TSP)，属于NP-hard。本文主要采取文献[5]中的贪婪算法进行求解，该算法属于组合优化问题的常规求解算法，虽然所得结果不一定为最优，但算法稳定可靠，具有实用性。具体实现过程见相关文献，此处不再赘述。

航线间转移时，转移路线与田块边界之间存在两种情况：一种是包含关系，即转移线段全部包含在某一田块的内部(包括边界上)，也就是无人机转移过程中始终处于边界内部；另一种是非包含关系，即转移线段上，除端点外，还有第3点处于边界外部，也就是说无人机转移过程中会跳到田块边界的外部。因田块边界外部可能存在高于作业高度的障碍物(尤其在田块间转移时)，为确保安全，有必要找出第2种情形下的线段并进行必要处理。

采用“安全边界相交性测试法”来判别某条候选转移路径线段Ji与田块边界多边形组{Bi}是否属于非包含关系。具体做法如下：根据障碍物情况，人为设立安全距离参数d和安全高度参数h；将田块边界多边形组{Bi}向外偏置d(外边界扩大、内边界缩小)，得到偏置后多边形组{B′i}；判断转移线段Ji与多边形组{B′i}的相交性，若相交，则视为“危险转移过程”。

如图6所示，对于田块F0的某条候选转移路径Ji(对应线段KiKj)来说，因为线段KiKj与外边界B0偏置后多边形B′0存在交点C1和C2，故该转移过程被视为危险转移过程。对于此类转移过程，实际操作时，飞机先从Ki点从作业高度h0拉升至安全高度h，转移至Kj点后再下降至作业高度继续作业。计算转移路线长度时在原计算值上加上2(h-h0)。上述过程中涉及到多边形偏置算法、多边形与线段间求交算法，也是计算机图形学中经典问题[21]，限于篇幅问题，不再赘述。

图6 安全边界相交性测试

1.7 航向优化

作业飞行的总线路长度和时间，除了与航线间调度次序有关，还与航线的方向有关。这里采用最小跨度法和步进旋转法进行综合优化，前者可精确确定单一凸多边形边界田块的最佳航向，后者则用于近似确定其他类型(凹多边形及孔洞多边形)田块的最佳航向。

区域跨度是指某一平面区域Q能够刚好被某一对方向为θ的平行线夹住时平行线间的间距D。此时，平行线称之为支撑平行线，θ称之为跨度方向角，D称之为区域Q在θ方向上的跨度。分析可知，凸多边形在任意方向上都有一对支撑平行线，且跨度值随方向角变化而变化，按照最小跨度对应方向布置填充线时数量最少，因此航线间转移距离将最短[24]。对于凸多边形而言，其跨度有3种形式，分别是点对式、点边式和边对式(图7)。文献[24]已证明凸多边形区域的最小跨度方向一定出现在点边式中。因此找出最小跨度方向角即为凸多边形田块的最优航向角。

图7 凸多边形跨度的3种形式

对非单一凸多边形田块而言，可采用“步进旋转法”进行近似求解。其本质上是一种穷举法，计算时，从0开始、以dθ为步距，计算所有可能航向下的航线，比较线路总长度，取最小值对应的航向θ为最优航向。

当θ≠0时，通过两次坐标变换可简化计算过程。先将多边形顶点绕其轴对齐包围盒(Axis-aligned bounding box，AABB)中心点(xr0,yr0)顺时针旋转θ(正变换)得到旋转后多边形顶点，基于活性边表法计算得到中间航线，再将航线端点(x0,y0)绕包围盒中心逆时针旋转θ(逆变换)即得到所需的真实航线的端点(x,y)。逆变换时计算公式为

(2)

而顺时针变换时，只需将θ值取相反数即可。

2 实例仿真计算

前述全部算法过程在Visual C++ 6.0平台上编码实现，在Intel Corei7-3370 CPU、3.4 GHz主频、8 GB内存、Windows 7操作系统环境下，分别对4组典型复杂性的人工假想边界和实际边界在多组不同工作参数条件下进行两次试验，以测试算法稳定性和计算效率。

实际田块1的边界利用Google Earth获取并以KML文件格式导出，实际田块2边界通过无人机开源软件Mission Planner获取并以其自定义的POLY格式导出，二者都经过坐标转换，将地球经纬度坐标转换成平面直角坐标后再进行计算测试[6]；假想田块边界利用CAD绘制并以DXF格式导入。计算过程中，算法对导入的边界自动进行合规性检查并进行分组；对于单一边界且为凸多边形的田块采用最小跨度法优化航向，其他的情形则采用“步进旋转法”(步距取0.5°)优化航向；采用贪婪算法进行航线排序优化，计算航线间转移距离时，若不考虑障碍物则直接计算直线距离，如考虑障碍物影响，则对每一潜在转移轨迹与安全边界进行安全检查，对不安全转移过程作高度方向上补偿，转移距离加上两倍高度差(即安全高度与作业高度之差)。结果数据中路径总长度通过轨迹线段长度求和得到，算法耗时通过GetTickCount()函数获取算法运行前后的系统时间后再求差得到，且测试时包括界面显示更新相关操作。

2.1 不考虑边界障碍物情况下的仿真试验

聚焦田块边界多边形的复杂性，不考虑边界上障碍物的影响或障碍物高度低于作业高度而不影响转移过程安全，测试结果相关数据如表2所示。表中，优化率ε为S2相对于S1的降低幅度，即ε=(S1-S2)/S1×100%，式中，S1为只进行航线排序优化而不进行航向优化(人工指定与边界上某一边的方向一致)时非工作路径(即航线间跳转路径)的长度，S2为同时进行航向及航线排序优化时非工作路径的长度。计算示例的部分结果截图(因篇幅限制，只取其中部分典型结果)如图8所示。

2.2 考虑边界障碍物情况下的仿真试验

为进一步扩大算法的应用范围，假设上述4个田块边界上存在障碍物且高度大于作业高度，在前述同时进行航向及航线排序优化试验基础上，继续增加转移过程的安全性判断及处理算法的测试。设置作业高度为2 m，安全高度为6 m，安全距离为1 m，测试结果如表3及图9所示。图9较图8增加了高度调整和安全边界，分别采用黑色小圆圈和粉色多边形表示。

3 讨论

不考虑障碍物影响时，从表2中数据及图8所示计算实例典型结果截图可以看出，本文算法能够处理各种复杂边界形状类型的田块，这些田块的边界形状复杂程度明显高于现有可查的相关文献中的算法实例，同时对于作业方向的优化是自动完成的，有别于现有商用无人机自带处理软件(包括Mission Planner)需要用户人工指定设置。

对比不进行航向优化而只进行航线排序优化的规划结果，二者同时优化时有效工作路径(即航线)的总长度基本不变，但对非工作路径的优化效果明显，优化幅度从23.04%到45.98%不等。这是因为有效工作路径要覆盖作业区域，其总长度与工作幅宽之积近似等于区域面积，而无效工作路径为航线间跳转轨迹，其总长度与航向和航线间调度次序高度相关。进一步，如不进行航线排序优化，则后果更加严重。如图8g所示，同样航向情况下，不排序情况下非工作路径长度为15 269.30 m，为表2中相同算例排序优化结果(对应图8e)2 058.88 m的7.4倍，是有效工作路径7 042.57 m的2.2倍。

表2 不考虑障碍物影响情况下的算法测试结果

图8 不考虑障碍物影响情况下部分计算实例结果截图

表3 考虑障碍物影响情况下的算法测试结果

图9 考虑障碍物影响情况下部分优化计算实例结果截图

此外算法耗时随着田块数量和边界顶点数目增加而增加，这是因为初步航线计算和航向优化都是以田块为单元逐一进行的，田块数量的增加必然会增加优化的工作量，而边界顶点数目的增加则会增加边界分组、多边形凹凸性判断及活性边表法求交处理的工作量。同时算法耗时还会随机器工作幅宽的减少和田块面积的增加而增加，这是因为幅宽越小、田块面积越大，航线的数量越多，需要求交的次数愈多、排序优化TSP的规模越大。对于航向的优化，基于最小跨度法的算法执行效率明显高于步进旋转法，但只适合于单一边界且为凸多边形的田块(如假想田块1)，通用性不如步进旋转法。算法总耗时从15 ms到19.2 s不等，考虑到计算优化过程是作业前离线完成的，对实时性要求不高，因此本文算法的实时性特性完全能够满足无人机使用时对路径规划的及时性需求。

考虑障碍物影响时，对比表3及表2中计算结果，在同样试验条件下，有效工作路径总长度几乎没有变化，而非有效作业路径长度普遍增加，主要原因是某些航线间转移路径中增加了高度调整的补偿而且航线间排序可能因此而改变。此外，算法耗时的增加幅度较大，特别是在多边形数目和顶点数较多且航线数量也较多时(如假想田块2在工作幅宽2 m时)，这是因为航线排序时要对每一潜在航线间转移过程进行安全性判断，极大增加了算法的工作量。好在算法的稳定性很好，计算效率也基本满足离线规划的要求。

无人机实际作业过程中，还可能遇到因为续航和载重等原因导致的消耗品按需补给、田块间高度差和多机协同作业等问题，这是路径规划工作中需要进一步深入研究的内容。

4 结论

(1)提出并实现一种对于田块边界形状具有普适性意义的旋翼无人机作业路径规划算法，可以快速获得凸多边形、凹多边形、带孔洞多边形或多个多边形等各种复杂边界田块下的飞行作业轨迹。

(2)利用假想田块、实际田块的多组仿真测试试验表明，所设计的算法运行稳定可靠、效率高、优化效果明显，可满足农用无人机作业时对路径规划过程的准确性、实时性和稳定性等相关要求。在不考虑障碍物影响时算法耗时15 ms～19.2 s，相比于只进行航线排序优化的情况，同时进行航向和航线排序优化后，航线间转移路径总长度下降了23.04%～45.98%，而考虑障碍物影响时的处理耗时也在离线应用的可接受范围内。