多约束情形下的农机全覆盖路径规划方法

陈 凯 解印山 李彦明 刘成良 莫锦秋

(上海交通大学机械与动力工程学院， 上海 200240)

0 引言

路径规划是现代智能农业机械实现自主导航移动和作业的三大关键技术之一[1]，决定农机工作总消耗并影响实际农机工作的作业效率和作业质量。

农机路径规划可以分为全局路径规划与局部路径规划[2]。其中全局路径规划是在作业环境已知的情况下，基于完全先验信息的路径规划，一般在农机作业前规划。局部路径规划是农机作业中采用传感器感知作业区域环境信息，以此进行的实时路径规划。全局路径规划常以减少作业总距离或时间消耗、提升作业覆盖率、提高农机作业效率[3]作为优化目标，优化方法可分为全局覆盖路径规划和全局点到点路径规划。全局覆盖路径规划主要以往复遍历平行作物行的方式进行作业[4-8]，通常在处理地块和生成作物行后，根据传统遍历方式确定轨迹。全局点到点路径规划是在已知作业区域环境信息后，规划从起点到目标点可行的无碰撞路径，目前比较成熟的算法有A*算法、蚁群优化算法、遗传算法、模拟退火算法和粒子群优化算法等[9-13]，这些算法可以根据不同优化目标进行起始点间最优路径搜索。

全局覆盖路径规划多针对地块边界约束和作物行约束，但大多使用传统的梭行、套行、开垄等[14-16]方式进行遍历，无法针对不同的农机转弯参数进行有效处理，对不同地块边界条件适用性较低。而全局点到点路径规划可以通过栅格化地块后，利用启发式算法等得到全局最优解，但仅限于多点间路径[17-19]，无法处理作物行约束，无法适用于需要往复遍历作物行的农业机械。

农业机械实现自主作业时需解决多种约束的影响[20]，处理好各项约束才能保证农机成功且高效地完成预定作业。常见的约束有农机地头转弯方式、地块边界、作物行行线、障碍物等。农机全覆盖路径规划可具体表述成基于以上的多约束条件下，以距离/时间/总能源代价为优化目标的路径规划。

农机作业总代价中作业代价与地头转弯代价占绝大部分，其余如出入地块等代价可忽略。当作业区域与作物行间距固定时，作业代价基本一致，因此该问题可近似为作业转弯代价的最优化问题。

本文提出一种考虑多约束的农机全覆盖路径规划算法，该算法将在量化不同地头转弯方式带来的代价，处理地块边界、障碍物等约束后，针对多种优化目标生成最优作物行。同时通过提出基于模拟退火和行数拆解与遍历合成的混合规则遍历顺序算法，实现大规模农机作业路径规划。

1 规划算法流程

本文提出的算法流程如图1所示，可分为多作业约束处理与最优遍历顺序求解两部分。作业约束主要包括农机动力学参数约束、农艺要求约束及地块边界条件约束，有效处理各约束可以保证作业轨迹按农艺要求正确完成作业。对于最优遍历顺序的求解，本文提出了基于模拟退火的混合规则求解方法，在大规模农田作业时对作业遍历顺序进行有效优化，降低各项作业代价。本文将二者进行结合，在处理约束的基础上优化了作业遍历顺序，实现了多约束下农机的全覆盖最优效率路径规划。

图1 全遍历作业路径规划流程图Fig.1 Process of coverage path planning

2 路径规划中的各种约束

2.1 地头转弯方式

农机的最小转弯半径与转弯间距之间的关系决定了其在地头换行时的转弯方式，一般可分为如图2所示的弓形、梨形、鱼尾形。图2中阴影部分为作业区域。当转弯半径为r，作业宽幅为w的农机在相邻作物行间换行时，由几何关系可知，当r小于或等于w/2时选择弓形转弯；当r大于w/2时选用梨形或鱼尾形转弯，两者可相互替代。

鱼尾形转弯虽然消耗小于梨形转弯方式，但过程包含倒车操作，对农机的转弯性能要求高，难度大，因此本文中路径算法优先考虑弓形、梨形转弯方式。

农机转弯的代价有：距离代价Cd，指单次转弯的总距离；时间代价Ct，指单次转弯的总时间；能量消耗代价Ce，指单次转弯的总能量消耗。农机转弯相对总作业时间的占比小，燃油消耗与时间近似成正比，因此可近似使用转弯时间代价表示。

图2 地头转弯方式Fig.2 Ground turning modes

实际作业时农机转弯性能受到不同农艺要求的约束。为实现农业作业无重复、无遗漏的目标，车辆在进行转弯时必须考虑到农机具所支持的转弯半径，如植保机等有较宽作业范围的机具需以较大半径转弯以避免重复喷药。因此实际转弯半径需综合考虑车辆动力学转弯半径rveh与机具转弯半径rag，本文取二者中较大值进行转弯代价计算。

根据各项农机参数，单次跨越n行的转弯距离代价Cd与时间代价Ct分别为

(1)

(2)

其中

r=max(rveh,rag)

式中vf——农机地头区域直线行驶速度

vt——农机转弯速度

α——当前作物行与地块边界夹角

车辆在不同作物行间的换行代价可使用转弯代价邻接矩阵L表示。对于作物行总数为N的地块，其邻接矩阵L大小为N×N，矩阵中元素Lij表示在第i行与第j行作物间进行单次转弯所消耗的代价，若i与j相等，则代价为0。某农机在某田地中5行作物行间的时间转弯代价邻接矩阵为

(3)

其最小转弯半径r与作业宽幅w的关系为w

2.2 地块边界

地块边界不仅是地块特征也是作业路径规划的约束，农机需在不越界的基础上完成作业。作业地块边界、地块内部的障碍物可由首尾相连的直线段构成的多边形来描述，一个完整地块可由一个地块外边界多边形和内部数个障碍物多边形组成。一般可使用卫星定位采集设备围绕农田地块边界或障碍物采样得到离散经纬度边界坐标点集，对采样点进行拟合可得到最终多边形边界。

本文基于道格拉斯-普克算法[21-22](Douglas-Peucker，DP)进行地块外边界拟合。DP算法可在保留原几何图像骨架的基础上进行线状要素抽稀实现多边形拟合，步骤如图3所示：在图3a中待处理点集的首末点间连一条直线，求其余各点与该直线的距离，找出最大距离dDPMAX及其对应的点P；设定阈值dthreshold，若dDPMAX

图3 道格拉斯-普克算法Fig.3 Douglas-Peucker algorithm

直接采用DP算法进行拟合时，若拟合后边界超出地块实际边界会导致农机越界。同理，若拟合后障碍物边界无法包络实际障碍物，会使农机与障碍物发生碰撞。本文针对边界拟合提出了内缩改进的DP算法，由于拟合后边界线与边界采样点距离最大为dthreshold，可在采用DP算法拟合后，将所有拟合后边界等距离内缩dthreshold，此时所有边界均在原采样点内部，确保所规划轨迹不越界。

对于障碍物边界，由于无需考虑形状影响，直接寻找障碍物采样点集的最小凸包即可完整外包络采样点。图4为对一含障碍物地块分别进行地块边界拟合与障碍物边界拟合后的结果，图中黑色边界为原采样边界，阴影部分为拟合后有效地块。

图4 边界内缩与外包络拟合Fig.4 Inner boundary contraction and outer envelope fitting

农机作业要求预留转向地块，该地块位于作物行两侧，用于农机完成转弯换行操作。如图5，一完整地块可分为工作区域W与转弯预留区域T，在其区域内分别进行作物行遍历与转弯换行操作。T宽度取决于转弯方式、作物行方向与农机自身参数，弓形转弯宽度dU和梨形转弯宽度dΩ分别为

dU=r(1+cosα)+w/2

(4)

(5)

图5 转向预留区域分割示意图Fig.5 Segmentation schematic of reserved turning area

已完成边界拟合的多边形地块可通过将各边向地块内部平行偏移边界的方式划分W与T区域。本文采用基于角平分的等距线法来实现该边界偏移：当各边均等距平移时，平移后两边交点位于原两边角平分线上；当各边平移距离不同时，可由不同边平行距离的比例来划分两边夹角，再求出平移后交点。

图6 角平分法偏移边界Fig.6 Parallel migration of angle bisectors

如图6所示，AB、AC两边界向内平行偏移距离分别为d1、d2，偏移后边界A′B′、A′C′交于点A′，其中AB、AC夹角为αA，有

β=αAd2/(d1+d2)

(6)

d′=d2/sinβ

(7)

式中β——角分线方位角

d′——角分线偏移距离

求解β与d′后可由点A坐标计算得到点A′坐标。对所有相邻边重复以上过程，得出所有工作地块边界顶点，按顺序连接后完成预留地块划分。

若实际作业条件允许，可选择进行封边操作，即在工作区域内正常作业完成后，遍历转向预留区域一周进行作业，进一步提高地块作业覆盖率。

2.3 非常规地块

农机作业需沿作物行作业且只能在地头区域转向换行，不能随意跨越作物行。对于常规凸多边形地块，MICHEL等[23]已经证明最佳作业方向应与多边形最长边平行，只需判断地块最长边，按该边角度划分等距平行线作为作物行行线即可。

相对一般凸多边形地块，非常规地块指含有障碍物或地块边界为凹多边形的地块。若使用上述寻找最长边的方法处理非常规地块则同一条作物行线可能会与障碍物或地块边界产生多个交点而带来额外的转弯次数，因此不能直接将最长边方向作为作物行倾角。

针对含障碍物地块，可将其分割为多个无障碍物子区域后，分别进行全覆盖规划后再连通各子区域，完成全地块的遍历[24]。一般可用梯形法[25]或线扫法[26]进行地块分割，本文使用基于最小高度和的线扫法完成区域分割。

2.3.1子区域分割

线扫分割法如图7所示，使用倾角θ的作物行从左至右划过整个区域，可按作物行与原地块、障碍物的交点情况，将地块分割成子区域①、②、③、④。此时全部子区域高度之和为

(8)

式中S(θ)——地块划分后总高度和

Hp(θ)——原地块相对θ角分割线的高度

Hi(θ)——障碍物i相对θ角分割线的高总高度和最小时该区域划分为最优划分。

图7 线扫分割法Fig.7 Line-sweep decomposition

2.3.2作物行线方向优化

当作物行间距固定时，转弯次数与子区域高度成正比，为此需找到最优作物行方向θ，以该角度线扫凹形地块或分割含障碍物地块，并按该角度在地块内划分作物行时，转弯次数最少，所生成作物行为最优作物行。

(9)

(10)

式中ns——地块及障碍物总边界数

θi——第i条边界角度

li——第i条边界长于θ′与θ相近，因此sin(|θ′-θi|)与sin(|θ-θi|)近似相等，且方向为θ时与对应平行边界无交点，显然此时N′raw>Nraw，即θ角对应作物行数为邻域内极小值。

图8 作物行数极小值示意图Fig.8 Sketch of minimum crop row number

2.3.3区域遍历顺序优化

完成地块分割与作物行线生成后，农机在分割后各子区块内的行走方式与无障碍物地块一致。在各子区域内，农机作业出口位置只由作物行数决定：偶数行数情况下将会在出发点同侧，奇数行情况下则相反。当车辆出入口有位置要求或农机在子区域的障碍物一侧结束工作时，受到作物行约束会难以驶出农田或前进到下一地块。可采用增加非作业路径的方法来调整出口位置。若某子区域存在n行作物行，定义遍历顺序集为

S={[s1,w1],[s2,w2],…,[si,wi],…,
[sn-1,wn-1],[sn,wn]}

(11)

式中si——农机路径中第i行在作物行中的顺序

wi——农机路径中第i行是否作业，若工作值为1，否则为0

当需调整出口位置时，在原n-1行空载不作业，在第n行正常作业后返回n-1行进行作业，此时出口在n-1行上，且位于原出口另一侧，即调整前后的遍历顺序集为

Sbefore={[s1,1],[s2,1],…,[sn-1,1],[sn,1]}

(12)

Safter={[s1,1],…,[sn-2,1],[sn-1,0],[sn,1],
[sn-1,1]}

(13)

以图9为例，原遍历顺序集为{[1，1]，[2，1]，[3，1]，[4，1]，[5，1]}，出口位于入口另一侧，当需要调整出口位置时，改变遍历顺序为{[1，1]，[2，1]，[3，1]，[4，0]，[5，1]，[4，1]}，此时通过增加非作业路径的方法改变总行数的奇偶性，使子区域出口在入口一侧。

图9 切换子区域出口位置示意图Fig.9 Sketch of changing subdomain exit position

增加非作业路径方法会额外带来单次作物行行走与转弯代价，但农机在非作业路径无需慢速作业，可快速通过节省时间代价。该方法以较小的代价解决了无法避免的车辆出入、切换地块问题，同时对作业条件和农机农具的行驶条件如倒车等无要求，具备一定的可行性。

3 基于多约束的遍历方式

现有农机作业行走路径模式常指梭行、套行、开垄形、闭垄形这些传统规则走法。下面进行最优作物行遍历顺序的生成研究，在上述各种约束条件下，求解以最小代价，无重复、不遗漏地遍历作物行的最优顺序。作物行遍历问题的本质为经过全部作物行且中间换行代价和最小的旅行商(Traveling salesman problem，TSP)问题[28]，其中换行代价使用转弯代价邻接矩阵求得。面对TSP问题，一些无规则启发式算法如模拟退火(Simulated annealing，SA)算法[29]理论上可以应对地块中的作物行约束。但当需面对大型农田时启发式算法面临规模效应的挑战，除了计算时间长外，其更容易陷入局部最优解，寻优率低[30]。为弥补标准模拟退火法的缺点，本文提出基于模拟退火优化的最小单元拆解与遍历合成的方法。

本文基于模拟退火的行数拆解与遍历合成的路径规划的实施路线与一般规则实施路径区别如图10所示。本文方法为根据当前作业参数的农机，通过改进模拟退火算法求解最优遍历单元。进行路径规划时将完整地块划分为数个重复最优单元，通过不断循环各单元，直至遍历完整块农田。因此，对全区域的遍历求解可转化为对最优行走单元的求解。

图10 路径规划实施路线Fig.10 Roadmap of path planning

3.1 基于改进模拟退火的最优单元求解

遍历最优单元为利用农机参数等先验信息在全覆盖规划前得到的可重复小规模农机路径。由于单元行数较少，避免规模效应的同时可充分利用模拟退火搜索的便利性完成求解。最优单元行数可通过求解最小平均单次转弯代价得到，对于单元内走法，本文采用引入寻优扰动算子的改进模拟退火算法进行求解。

模拟退火算法通过模拟金属退火的过程，从某一较高初始温度To开始降温，不断改变当前解s，若其评价值F(s)变小，则接收新解，否则以概率PSA接收新解，当温度降至终止温度Tf，结束降温。多次退火后所获稳定解即为最终寻优解。对于本文所述问题，所求解s为单元内遍历顺序，si为单元内第i行对应实际作物行顺序，评价值F(s)为单元内以及前进到下一单元的总转弯代价，其定义为

(14)

式中Nu——单元内行数

Lij——转弯代价邻接矩阵中i行j列对应元素值

为改善算法性能，本文引入了寻优扰动算子，在产生新解时根据当前温度高低调整当前需调整路径行数，进一步提高在高温时算法的搜索范围而不影响低温时的收敛寻优。每次调整的路径行数Nadjust为

(15)

式中 int()——取整算子

Nmaxa、Nmina——允许最大和最小的调整行数

T——当前温度

由此，使用该改进模拟退火算法求解最小单元走法步骤为：

(1)设定初始温度To，根据当前单元行数产生随机初始解s，计算其转弯代价评价值F(s)。

(2)降低温度，根据当前温度T计算寻优扰动算子，随机选取Nadjust行，对其遍历顺序进行扰动交换，产生新解s′并计算当前新的评价值F(s′)。

(3)新解s′的接受：若F(s′)

(16)

(4)判断算法终止条件，若当前温度T

模拟退火算法所求得的无规则走法，没有改变车辆行驶条件与作业条件，仅对车辆在转弯区域内的换行顺序进行优化，不改变全覆盖作业过程，对农艺无影响，具备可行性。

3.2 基于最小单元的行数拆解及遍历合成

在单地块内，由于无法事先确定作物行数，当总行数不是最小单元行数的整数倍时，单元重复遍历后会存在多余行。因此本文通过行数拆解与遍历合成的方式完成规划，拆解前后总行数与单元行数的关系为

(17)

其中

k=int(Ntotal/Nunit)

式中Ntotal——地块总行数

Nunit——最小单元行数

Nleft——单元遍历后的剩余行

k——完整单元数

根据式(17)，任意总作物行数均可拆解为数个单元行数Nunit与剩余行数Nleft，只需分别对各单元行及剩余行内的走法进行规划即可完成全地块遍历。由此本文提出最优路径集的概念，最优路径集包含从第1行至2Nunit-1行各情况下的最优规划走法，通过改进模拟退火算法求解各情况下走法。最优路径集在全覆盖规划前生成，全覆盖规划时可按式(17)对总作业行数进行拆解，再于最优路径集中分别提取单元内走法与剩余行内走法，通过合成即可生成完整地块轨迹。

3.3 行数拆解及遍历合成应用

根据该行数拆解及遍历合成方式，本文进行了不同农机参数条件下的最优单元和路径集生成。本文采用±Ui、±Ωi、±∝i分别表示隔i行的弓形、梨形、鱼尾形转弯，其中±表示正反向。采用F(Ui)、F(Ωi)、F(∝i)表示对应的单次转弯代价，可由F(s′)

(1)r

此情况下转弯代价关系式为F(U1)

图11 梭行遍历Fig.11 Fusiform traverse

(2)w/2≤r

此情况下转弯代价关系式为F(U2)

(3)w≤r<3w/2(F(Ω2)

此情况下转弯代价关系式为F(U3)

(4)w≤r<3w/2(F(U4)

此种情况下的转弯代价关系式为F(U3)

(5)3w/2≤r<2w

此种情况下的转弯代价关系式为F(U4)

图12 最优路径集Fig.12 Aggregation of optimal path

以上已涵盖大多数农机与作物行参数关系，因此其余情况本文不再展开。

4 测试与分析

4.1 测试方式

在北京小汤山国家精准农业研究基地对该规划算法进行了验证。地块数据使用卓林A5型多功能测亩仪的点采集功能绕基地中形状不一的多个地块各一周采样得到，平均水平定位精度为0.55 m。测试农机设备为约翰迪尔6B-1204型拖拉机，其最小转向半径为5.2 m，拖挂耕地设备作业宽幅为6 m，工作行进速度为2 m/s，转弯平均速度为1 m/s。配置编程环境为Windows 10、Visual Studio 2019。

4.2 评价指标

农机路径可使用总时间代价CT、总距离代价CD、作业覆盖率Ccov、工作占空比D等指标来评价其质量。其中，总时间代价CT与总距离代价CD分别为农机沿规划路径完成遍历所需的总时间与总距离，表述农机完成工作所需的总资源，代价越低，该轨迹质量越高。

工作占空比D表征作业时间与总时间的比例，工作占空比越大，农机工作效率越高，公式为

(18)

式中twork——农机负载工作时间

tall——总作业时间

作业覆盖率Ccov表征有效作业面积在地块总面积中的比例，覆盖率越大，土地利用率越高，规划质量越高，公式为

(19)

式中Lw——工作路径总长度

Sfield——地块有效总面积

4.3 结果与分析

为验证本文提出的多约束处理方案的正确性，按不同约束条件在北京小汤山国家精准农业研究基地采集了多个地块边界数据，包括规则凸多边形地块(地块a)、不规则凹多边形地块(地块b)及含障碍物地块(地块c)。使用上述农机，按本文所述全覆盖路径规划方法对不同地块分别进行了规划，生成路径结果如图13，并求出各路径的作业覆盖率及作业占空比。考虑到农业生产可在作业完成后通过封边作业对转弯预留地块进行遍历以提高土地利用率，因此同时计算出封边操作后的作业覆盖率。结果如表1所示。

图13 地块边界及规划路径Fig.13 Plot boundaries and planned paths

表1 不同地块所得轨迹评价指标Tab.1 Trajectory evaluation indexes obtained from different plots %

由图13可看出，各路径均预留了边界转弯地块并根据地块形状生成了最优作物行，对含障碍物地块进行了区域分割处理，并针对多种转弯方式约束利用本文单元拆解合成方法，优先选取弓形转弯生成了最优单元，完成了全覆盖遍历。

分析表1可知，本算法面对包含不同约束的地块，作业参数均可满足规划需求，封边作业后各测试地块土地利用率均可达98%以上，最高可达99.87%。同时，作业占空比均较高，保证了农机的作业效率。说明本文提出的规划算法可以在有效处理各种约束的基础上完成农机全覆盖路径规划。

为验证本文提出的行数拆解及合成的遍历算法的性能，对规则凸多边形地块a分别使用传统规则遍历方法、经典模拟退火算法(取To=200℃，Tf=50℃，经200次重复退火取代价最小值)与本文提出的遍历方法进行规划，求出各方法所得轨迹评价指标，包括总时间代价、总距离代价、作业覆盖率(均不进行封边)、工作占空比。为验证不同农机参数对算法的影响，更换农机拖挂作业设备，改变作业宽幅为3.5 m，农机转向参数不变，再次规划并计算评价指标，结果如表2所示。

为验证不同地块约束对算法的影响，在同样的作业参数下使用各遍历方式对地块c进行规划。相对地块a，地块c更符合大田的特征：面积更大，作物行相对地块宽度更长，作物行约束更明显，同时存在障碍物约束，因此可用于检验算法对不同地块的适应性。地块c条件下各算法所得轨迹评价指标如表3所示。

表2 地块a不同走法下轨迹评价指标Tab.2 Trajectory evaluation index of plot a under different moves

表3 地块c不同走法下轨迹评价指标Tab.3 Trajectory evaluation index of plot c under different moves

分析表2可知，传统规则走法中：套行走法与闭垄走法均只适用于农机转向性能较差及地块宽度较小的情况，因此梭行走法各代价均相对较小，但逐行梭形遍历无法避免使用梨形转弯，需要预留更宽的转弯区域，作业覆盖率较低。非规则走法中：模拟退火算法结果与梭行法相近，但当作物行数增多时，规模效应使得模拟退火算法难以逼近全局最优解。而本文提出方法在模拟退火的基础上，降低了规模效应，使得结果更加接近最优解，相比传统规则走法节约距离消耗最高达30.3%，相比经典模拟退火算法，节约距离消耗6.9%，同时提高了作业覆盖率与作业效率。

对比分析表2与表3，作物行长度与地块宽度的比例增大，使得地块c相比地块a中各算法所得轨迹作业占空比均有提升，但传统规则算法受边界条件变化影响更大，鲁棒性低，未能有效处理不同地块边界约束。而本文提出的行数拆解及合成的遍历算法在不同地块边界条件下均可保持较高的作业占空比，在地块c中高达94.72%。说明本文提出的遍历顺序算法可以有效克服各约束的影响，以较优的代价和作业效率完成农机的路径规划工作。

5 结论

(1)针对大型农机作业路径规划中存在的多种约束，本文通过引入转弯代价邻接矩阵量化了不同地头转弯方式对规划的影响，采用内缩改进的道格拉斯-普克(Douglas-Peucker)拟合算法和采样点的最小凸包法分别处理地块边界和障碍物边界，避免车辆越界和与障碍物发生碰撞。使用角平分线的平行偏移法求得转向预留地块后，以转弯代价最小为优化条件对多种形状地块进行了最优作物行生成。

(2)在作物行的最优遍历顺序求解的问题上，提出了基于模拟退火的行数拆解及遍历合成的混合规则优化算法，该方法适应性强且优化高效稳定，解决了传统无规则模拟退火算法因规模效应而容易陷入局部最优解的问题。同时相对传统规则遍历方式，有效处理了不同边界约束，大幅减少了作业消耗，提升了作业效率。

(3)通过实地测试，测量了多个地块参数并使用本文方法进行了规划。所得规划轨迹的各项指标中，总距离、时间代价均低于传统规则走法与经典模拟退火算法，较传统走法最多可节约距离消耗30.3%，较模拟退火算法节约距离消耗6.9%。本文算法所得轨迹平均作业覆盖率达90.78%，平均作业占空比达85.10%，封边作业后地块有效覆盖率达98%，可以满足实际规划需求。