肖菲菲, 蒋蘋, 胡文武, 廖荣华, 张丹慧, 金生

(湖南农业大学机电工程学院, 长沙 410128)

目前，精准农业是农业智能化发展的研究热点之一，其中田间自动导航技术是完成精准农业的关键环节[1]。农业机械的自动导航技术现已被广泛应用于喷药、施肥和除草等农业作业过程[2-3]，该项技术不仅可以降低劳动强度和农业生产成本，还可提高农机作业效率。而在农业机械的自动导航过程中，农作物的苗带识别和定位是核心技术之一，在很大程度上影响农业机械的导航效果和作业的精准度[4-7]。苗带是指根据农作物种植方式，在同一方向上农作物连续性的集合，在作物生长初期，苗带为条状；在作物生长后期，作物封行后，苗带为块状。苗带识别是指通过某种方法区分农作物不同生长状态下的苗带区域与相邻苗带间的行间间隙所在位置。通过苗带识别，在可行走区域规划行走路径，使机械行走过程中避免碾压苗带，实现田间管理机械自主低伤苗作业，对促进精准农业发展有重要意义[8]。

当农业机械行驶在环境复杂的农田中，为了避免障碍物和压伤农作物，需要通过传感器进行感知识别，当前常用于苗带检测识别的传感器有视觉传感器和激光雷达等[9-14]。王姗姗等[15]通过视觉传感器检测水稻秧苗，采用基于特征点邻域 Hough变换算法完成水稻苗中心线提取，获得平均识别准确率为92%，其识别率较好，但响应速度慢，只适用于速度低于1.5 m·s-1。张博立等[16]利用机器视觉提取玉米苗带中心线，采用动态网格与分区域聚类算法相结合进行识别，获得平均玉米苗带识别准确率为93.4%，有效地处理了杂草对识别效果的影响，但由于图像相关处理问题，导致中心线检测失败,且目前仅适用于识别生长初期玉米作物行。姜国权等[17]通过粒子群优化算法获取麦田作物行的特征点，利用最小二乘法进行线性拟合，运行速度相对传统Hough算法有很大提高，但并不适用于处理杂草较多的图像。Hu等[18]利用二维伺服平台，通过激光雷达扫描获得水稻苗带点云信息，设计基于斜率虚化均值迭代聚类算法，采用斜率虚化处理点云数据，通过区域迭代聚类获得苗带特征，利用最小二乘法进行苗带轨迹线性化拟合，经试验获得横向最大误差为45 mm，中位点偏差为7 mm，但未考虑激光扫描偏移问题，且存在重复运行算法时出现聚类失败的情况。

由以上研究可知，视觉传感器成本较低，可实时获取作物行相对位置，因此广泛用于作物苗带识别，但传感器易受外界环境因素的影响，且数据处理相对繁杂，部分系统并不能满足真实的农田环境。而激光雷达具有易操作、分辨率高、抗干扰性强等特点，更宜用于农田环境，但目前大多研究集中于室内环境或已知结构的室外环境，对于水田的作物识别研究相对较少，且并未考虑扫描偏移对识别结果的影响。

综上，本文提出了基于偏移补偿模型的极值点聚类苗带识别算法，重点研究激光扫描偏移对识别结果的影响。首次提出了苗带识别的相关定义，根据获得的苗带点云线特征，采用极值检测算法获取聚类要素点，根据苗带特征，分析激光测距仪的扫描特点，建立偏移补偿模型，对获取的聚类要素点进行偏移矫正，最后通过最小二乘法进行线性拟合，获得苗带轨迹信息。提高了轨迹重合精度，增强了算法运行效率，获得的苗带轨迹参数可为田间管理机械低伤苗直线行驶提供理论依据。

1 材料与方法

1.1 系统构成及工作原理

导航定位控制系统由植保机行驶控制单元(electronic control unit)和苗带识别系统构成，其中苗带识别系统包含激光雷达、舵机调平单元与惯性测量模块，如图1所示。

注：1—施肥斗; 2—发动机; 3—苗带识别系统； 4—植保机行车控制单元；5—药箱；6—喷杆。

激光雷达安装于舵机平台上，惯性测量模块采集姿态数据，舵机调平单元通过采集惯性测量模块的数据，控制舵机进行俯仰、翻滚、旋转动作，进而实现激光雷达底座的俯仰角、偏航角、横滚角可控调节，并通过CAN接口与植保机ECU单元进行数据交换。ECU单元采集激光雷达扫描的点云数据与舵机平台角度数据后，通过处理获得水稻苗带信息，并识别可行驶区间，控制田间管理机械行驶于苗带行间间隙中，减少作业过程中对水稻的损害。

1.2 苗带识别系统硬件与数据采集软件

水稻苗带识别系统的基本结构(图2)由三部分构成，分别为激光雷达、二维舵机调平单元和支撑平台。本研究选取日本HOKUYO的UTM-30LX二维激光雷达，测量范围为100～30 000 mm，100～10 000 mm的测距误差为±30 mm，扫描角度范围-135°～135°，扫描周期40 Hz，角度分辨率为0.25°。舵机调平单元采用达盛舵机RDS3115，实现X、Y轴的水平调节控制，在履带底盘底部安装旋转调节舵机，实现X、Y轴的旋转调节，完成激光雷达俯仰角、横滚角、航向角的调节控制。为实现舵机角度初始化校准和激光雷达中位调整，采用核心处理器STM32F103配合MPU6050六轴姿态传感器，实现舵机角度调节。

注：1—控制单元; 2—支撑平台; 3—旋转调节舵机； 4—激光雷达；5—俯仰调节舵机；6—翻滚调节舵机。

激光雷达通过软件UrgBenri Standard V1.8.1进行数据采集，数据存储格式为ubh，使用Lenovo CPU1.60GHz在Win10环境下运用Matlab R2016b进行数据分析与算法研究。

1.3 点云数据获取与苗带轨迹提取分析

1.3.1点云数据获取为确定激光束发出位置，确保激光雷达采集数据的准确性，本研究在数据采集前通过最小二乘平差法求得真实距离与激光雷达所测距离之间的数学关系[19-20]对测量仪进行标定，以消除激光雷达自身扫描误差。苗带数据采集于2020年7月湖南农科院试验田中，此时水稻处于移栽后15 d，当天天气预报风速为2 km·h-1，环境温度32 ℃，目测水稻叶片处于静止状态。采用卷尺测量水稻高度、行株距和激光雷达距离水田平面的高度等参数，采用5组连续15株水稻苗信息求均值的方法，获得水稻高度和行距株距。通过测试计算，得到水稻苗与水面高度平均值为360 mm，平均行距为220 mm，激光雷达扫描中心点到水田水平面的垂直距离为960 mm。在实际稻田环境数据采集过程中，二维舵机调平单元的俯仰扫描初始角度为9°，步进角度0.5°，终止角度17°，激光雷达扫描角度范围为-16°～16°，扫描宽度对应5行苗带，以此试验条件获取苗带点云数据。

1.3.2坐标转化激光雷达测距显示参数主要为当前激光雷达扫描的起止角、步长、扫描线数，导出数据为距离与角度的极坐标数据，系统参数还包括当前激光雷达的俯仰角度数据。图3为苗带识别系统工作时的相关参数，以水田水平面为参考，以激光雷达激光旋转中心为原点，垂直方向为Z轴，激光雷达初始光轴与机具直行重合方向为Y轴，垂直前进方向为X轴，激光雷达俯仰角度θ，激光测距长度为L，对应激光线扫描线角度为γ，Y轴对应激光雷达中心位置偏移初始角度为γ0。为方便数据分析，对激光极坐标数据进行三维坐标转换，基于上述坐标转换模型获得被测点相对激光雷达当前位置的相对直角坐标表达式(1)。

图3 扫描系统相关参数

(1)

1.3.3苗带轨迹获取方法为获得最终苗带轨迹，本研究系统将激光雷达和惯性测量模块采集的数据融合，获得被测目标点的位置坐标，运用Matlab R2016b进行数据分析与算法研究。根据农业机械作业环境和激光雷达点云数据特征，确定本文数据处理流程为数据预处理、点云极值点提取、苗带特征点聚类和轨迹获取与评价四个部分。数据预处理主要包括点云数据格式转换、数据去噪和数据坐标转换等步骤；点云极值点提取主要实现苗带特征点的提取；偏移补偿均值聚类主要进行苗带特征点聚类；轨迹获取与评价主要完成聚类点的线性拟合，并对获取的轨迹线进行有效性评价。

1.4 基于偏移补偿模型的极值点聚类

1.4.1苗带模型建立由于田间管理机械直线行驶仅需要获取以水平面为参考的平面坐标，因此，对三维数据平面化，获得基于激光雷达为原点的XY坐标定位数据。

(2)

为了方便计算苗带位置，获得苗带函数，根据激光雷达扫描测距方法建立理想苗带模型，如图4所示，长条框为苗带，激光测距仪距离植株距离b，苗带宽度d1,相邻苗带间距为d2。为验证模型的正确性，使用5个纸箱进行模拟试验，中心间距500 mm，UrgBenri Standard V1.8.1采集图像。由图像对比可知 ，建立的模型与模拟实验获得的点云线相对应，图像在X轴方向都具有周期性，且为具有规律性的波浪峰，模型结果成立。

图4 数据采集模型

从数据采集模型(图4)可知，红色线为激光雷达采集的点云数据线，图像在X方向具有周期性，周期为d1+d2，且图像关于Y轴对称，根据模型图和已知参数计算模型函数，函数共分为三段，计算得出函数关系式(3)。

(3)

1.4.2建立偏移补偿模型根据建立的模型可知，由于苗带的物理遮挡，导致部分位置的点云数据无法正确采集。在扫描中间苗带时(图5)，扫描线经过Q点后，右侧面的点云数据由于遮挡无法采集，使得在QN线上的点云数据产生了偏移，所以需要建立偏移补偿模型进行矫正。在激光扫描范围内，以原点为中心左右两侧最高峰为大波峰，最低峰称为小波峰，Y小于小波峰以下的数据相对较为准确，偏移量较小，可通过聚类算法直接获得苗带聚类点；Y在小波峰和大波峰之间的数据偏移较为严重，需对获得的苗带聚类点进行偏移矫正。基于以上问题，建立偏移补偿模型，如图5所示，在纵向小波峰到大波峰范围(GA≤Y≤NA)内，获得偏移量ΔX的矫正范围。

图5 偏移补偿模型

在EM线上，通过苗带模型函数，在点E、M处分别得到如下公式。

(4)

(5)

(6)

在QN线上，通过苗带模型函数，在点H、N处分别得到如下公式。

t×tanγ1=GA=MB

(7)

(8)

(9)

则苗带间隙最大偏移量计算如下。

(10)

因此，偏移量范围为0≤ΔX,基于此建立的偏移模型在同一波中，取低于小波峰的点云偏移量趋向0，在小波峰与大波峰之间的点云偏移量为ΔX，获得矫正模型(图6)，表达如式(11)所示。

图6 矫正模型

(11)

1.4.3聚类要素点的获取根据建立的模型(图4)可知，每条点云数据线在X方向都存在规律性的波浪峰，而波峰波谷所在位置与苗带和苗带行间间隙相关，然而波峰波谷位置又对应函数中的极值所在位置，因此，可通过计算函数的极值所在位置得到与苗带位置相关信息。本研究采用极大值检测算法进行极值检测，寻找苗带位置，由式(12)确定。

s=find{diff[sign(diff(w))]<0}+1

(12)

式中，diff指差分运算，sign为符号函数，find表示查找点云数据线中满足条件的数据位置[21]。

当采样点满足式(12)时，其点云数据在邻域范围内呈先增加后减小的趋势，即局部极大值点，获得波峰位置；将原数据的纵坐标值取反，再利用函数查找波谷位置。因农田环境复杂，水稻叶片具有下垂交叠等特征，且此特征也可作为判别苗带位置相关因素，为避免信息缺失，所以阈值范围取最小值到最大值，使用Matlab 绘制多角度点云数据线的波峰波谷图，进而获得苗带聚类特征信息。

1.4.4改进K-means均值聚类均值聚类算法的基本思想是以空间中k个点为中心聚类，对最接近他们的对象进行归类，通过迭代的方法不断调整各分类中心位置，直到收敛[22]。聚类算法中聚类数量与初始聚类中心的选择会直接影响聚类结果，通过改进K-means聚类算法提高苗带点云数据的聚类效果。通过矫正函数获得最佳聚类个数K,具体流程如图7。初始中心的选择通过密度公式(13)确定，在r∈(0,d1+d2)范围内，定义对应运算区域的聚类点位置X，计算xi内的点密度ρ(xi),当密度达到最大时，在对应范围内采用欧式距离度量方式选取距离xi最远的点作为第(i+1)类初始聚类中心xi+1，并且除去xi、xi+1，以此类推，重复上述步骤再获得下一个聚类中心，不断迭代，直至获得第K个聚类中心结束。

图7 获取聚类点数流程

ρ(xi)={x∈X|dist(xi,x)≤r}

(13)

1.4.5苗带轨迹处理为便于获取苗带轨迹，获得更多用于线性拟合的聚类点，同时降低数据采样维度，采用纵向等区间划分增加数据获取量，获得更为准确的轨迹参数，对Y轴向范围进行n等分区间聚类。通过分层聚类后，获得分段范围内的聚类点。

(14)

式中，k为单次聚类数量，(x,y)为聚类点对应的坐标位置,n为Y轴等分处理的次数，从而获得nk个聚类点。则对应苗带行间间隙点为(xv,yv),矩阵为O(n)。

(15)

(16)

(17)

(18)

基于偏移模型进行矫正获得聚类矫正后的矩阵T。

(19)

从而获得对应Y轴向的s组聚类点矩阵M(s)。

(20)

对上述矩阵坐标点进行线性拟合获取对应苗带行间间隙的轨迹方程。对聚类点坐标矩阵通过最小二乘法进行线性拟合，获得Ak和Bk的值，从而获取苗带对应的理论轨迹。

(21)

(22)

对应苗带间隙轨迹表达式如下。

y=Akx+Bk

(23)

1.5 苗带轨迹评价

参考文献[23]的评价指标，在获得理论苗带轨迹后，与对应水稻苗带真实轨迹进行对比，采用测量轨迹与真实轨迹平行度和中心偏差距离两个指标评价算法的平行度(式24)和中心偏差(式25)。

Err_max=|xk-x′k|max,yk∈(a～b)

(24)

Err_mid=|xk-x′k|max,yk=(a+b)/2

(25)

式中，xk为对应第k行x的真实值，x′k为对应第k行x的测量值。

2 结果与分析

2.1 测量仪标定误差分析

通过5组不同距离的数据，对激光测距仪进行标定，标定后实际距离与激光雷达测量结果见表1。可以看出，随着激光雷达所测距离由0.645增加到14.389 m,数据标定前与标定后的误差也逐渐增加，但标定后的误差明显低于标定前的。当被测物体距激光雷达14.389 m时，标定后的误差为0.011 m，相对标定前的误差0.037 m明显减少。说明测距仪使用前有必要进行标定。同样，完成二维舵机平台与激光雷达的联合标定。

表1 雷达标定前后测量值误差

2.2 试验数据模型构建

农田环境中，作物通常采用彼此平行的种植方式[24]，采集试验区域对应有效测量范围的水平投影数据。以激光雷达为坐标原点，为消除侧偏导致的方位差，激光雷达前视方向与中间列位置基本平行，为满足前视行驶需要，设定前视水平投影距离为3 030 mm，前视水平投影探测有效范围为2 000～3 030 mm，有效长度为1 030 mm，左右水平投影探测有效宽度为-500～+500 mm。测试环境的真实苗带轨迹如式(26)。

(26)

则测试环境的可行驶苗带行间间隙轨迹如式(27)。

(27)

2.3 苗带点云数据分析

2.3.1苗带三维点云显示及分析通过Matlab R2016b对激光测距仪采集的点云数据进行预处理，绘制苗带轨迹三维立体图，如图8所示。

图8 MATLAB转换后的点云数据

根据田间采集的实况与三维点云图外观对比分析可知，对应的苗带信息清晰有规律，且相邻苗带行间间隙的初始轨迹鲜明。箭头标明了苗带行间间隙轨迹所在位置，但受叶片倾斜折弯和交叉下垂以及激光雷达扫描方式的影响，行间间隙宽度较窄，且随着前视距离的增大间隙轨迹不易识别。多角度扫描轨迹线在靠近X轴有五个突出的波谷，为激光雷达扫描的第一行水稻苗所呈现出的形状，且每根点云数据线都存在这种规律性的波浪峰状特征。

2.3.2点云极值分析通过对每根点云线进行极值检测，获得波峰波谷数据点位置，获取聚类要素点，如图9所示。

图9 波峰波谷数据点

注：图中不同颜色的圈圈代表不同的分类，*代表此类型的聚类点。

表2 最终聚类中心

2.3.4苗带轨迹提取通过聚类后进而获得水平矩阵R(n)。

(28)

(29)

则对应苗带行间间隙点为(xv,yv),矩阵为O(n)。

(30)

(31)

(32)

对应获得水平聚类点矩阵Q。

根据实际采集情况，基于偏移模型计算获得偏移量范围为0≤ΔX≤174 mm，获得聚类矫正后的矩阵T。

(33)

通过平面投影处理，获得如图11的分段均值聚类分布图，其中黑点为聚类中心点，红点为矫正点。

注：图中不同颜色的圈圈代表不同的分类，*代表此类型的聚类点。

获取Y轴向对应的4组聚类点矩阵M1～M4。

(34)

对应上述纵向4组苗带行间间隙聚类点，在Matlab下进行线性拟合处理，获得轨迹矩阵如式(35)，对应轨迹如图12所示。

图12 苗带行间间隙轨迹对比

(35)

2.4 轨迹结果分析

在y∈(2 000,3 000) 范围内，计算获得轨迹与实际轨对应横向最大偏差Err_max和中位点偏差Err_mid。

Err_max=[16111614]

(36)

Err_mid=[6461]

(37)

表4分别统计了采用两种不同算法对苗带点云数据处理后的X轴最大横向偏差、中位点最大偏差和平均耗时。

表4 算法检测偏差及检测时间对比

从表5中数据可得出，在前视距离2 000～3 000 mm范围内,通过斜率虚化聚类算法对点云数据进行处理，获得轨迹与实际苗带轨迹X轴横向最大偏差为45 mm, 中位点最大偏差7 mm，平均耗时1.64 s。本文采用基于偏移补偿模型的极值点聚类苗带识别算法获得X轴横向最大偏差控制在16 mm以内，中位点最大偏差为6 mm，运算耗时1.50 s，相比斜率虚化聚类算法X轴向最大偏减小了29 mm，中心点最大偏差减小1 mm，耗时减少0.14 s，检测效果明显提高，所检测苗带行间间隙轨迹与实际轨迹具有较好的重合度，满足水田管理机械免压苗行驶的需求(轮宽10 cm)。

3 讨论

苗带轨迹的识别是实现田间自动导航技术关键。现有对作物的检测大多针对旱地作物,如蔬菜、玉米、小麦等，对于水田的作物识别研究相对较少。本文首次提出了苗带识别的概念，利用二维激光雷达与舵机调平单元相结合，以获得三维扫描功能完成苗带点云数据获取，进行水稻苗带识别研究，避免了视觉传感器易受外界环境影响，数据处理的算法复杂度较高等问题，重点研究激光扫描偏移对识别结果的影响，为解决激光雷达扫描偏移和提高轨迹识别精度等问题，提出了基于偏移补偿模型的极值点聚类苗带识别算法。根据获得的点云线特征，利用极值检测算法获得聚类要素点，分析激光雷达的苗带扫描特征，建立了偏移补偿模型，结合改进的K-means聚类算法进行偏移矫正聚类，最后通过最小二乘进行线性拟合，从而实现苗带轨迹可行驶区域的识别。

通过对前视距离为2 000～3 000 mm范围内的苗带轨迹进行识别，获得最大横向偏差低于16 mm,中位值偏差控制在6 mm以内。相对斜率虚化聚类算法，考虑了激光扫描偏移等问题，为提高聚类效果，改善原始K-means聚类算法中聚类数量与初始聚类中心需自行设定等问题，结合苗带特性对其进行了改进，试验结果表明，基于偏移补偿模型的极值点聚类苗带识别算法轨迹识别精度更高，平均运算耗时更少，算法鲁棒性更强，可为田间管理机械进行低伤苗自主导航行驶提供理论依据。但本文算法进行了Y轴纵向等分聚类，在选取不同等分数量时，可能会出现聚类效果不佳现象，后期需进一步改进。文中采集数据的试验田为机插秧的标准化水田，测量对象无缺失，测量现场环境良好，但实际农田环境中有诸多不确定性因素影响，需要在后期进行进一步的研究。