近十年来，我国水果产量、消费数量与种植面积均居世界首位^[1]，从品类来看，苹果产量占比位居我国水果产量总占比前列，且其种植面积与产量呈持续增长趋势。果园种植已成为我国继粮食、蔬菜之后的第三大农业种植产业，而生产模式和生产成本制约着水果产业的发展，优化生产模式、提高生产效率^[2-3]和节约生产成本离不开农业机械装备的智能化与无人化。

可行驶区域精准识别是农业机械装备实现智能化的关键技术之一，可为后续采摘、喷药以及其他果园管理工作提供可靠基础。现有研究多以激光雷达和视觉传感器识别为主，激光雷达适合全天作业且抗光照干扰程度较强，但其价格高昂、易受气候影响，且工作时被测物体与雷达间不能有障碍物^[4]；视觉传感器灵活性大、获取信息完整丰富且成本较低，符合农业装备经济性要求^[5]，但其易受复杂环境的影响，如光照、阴影以及果树行间杂草等^[6]，机器视觉发展潜力较大，正逐渐成为可行驶区域识别与路径导航的研究热点^[7]。基于机器视觉识别可行驶区域在果园中的应用可分为以天空为背景和以地面为背景两类，Radcliffe 等^[8]提出以果树树冠间天空为背景，使用绿色分量(G分量)增大树冠与天空的对比度进行路径平面的提取，通过天空质心计算车辆的航向，但该方法要求树冠饱满且为绿色。目前，国内已有研究多以地面背景为主，如林桂潮等^[9]根据道路特征提出果园道路边缘提取(Road edge extraction，REE)算法，利用梯度幅值法进行阈值分割并结合灰度值对比度约束区分道路左右边界，改进最小二乘法和随机采样一致性(Random sample consensus，RANSAC)算法与给定多项式模型相结合得到拟合后的道路边界；聂森等^[10]利用HSV颜色空间中H-V灰度化因子对图像进行增强，并使用Hough变换对果树底部靠近地面处树冠进行拟合得到树行直线，但该方法在光照较弱以及树叶凋落季节效果较差； 通过提取树干与地面的交点作为特征点，形成可行驶区域的左右边界线，但当杂草生长较高并遮挡树干与地面交界处时，无法精确提取特征点，且果树阴影也会对特征点提取产生影响^[11-12]。综上所述，目前已有的以天空为背景进行区域识别与路径提取的研究多对树冠要求较高，当果树处于凋落期枝叶稀疏且颜色特征并不明显时，特征提取以及拟合效果较差；以地面为背景时提取效果易受光照、树冠阴影以及杂草等因素的影响。

本文提出一种复杂环境下以天空为背景的果树行间可行驶区域识别算法，通过天空与树冠交界处直线获得可行驶区域边缘直线斜率信息，利用剔除后特征点信息和剔除阈值获得关键点坐标，以斜率为约束条件，代入关键点得到可行驶区域边缘直线方程并使用最小二乘法进行拟合。

1.   树行预处理

本文从山东理工大学与山东中以现代智慧农业有限公司共建的标准化苹果园实际环境出发，果树行距约4 m，果树高度3.7 m左右，树行长80 m。根据果实生长周期对复杂环境进行定义，即由果实生长期内地面上掉落的果实、树冠阴影，非生长期内地面落叶、树干阴影、树干间背景以及全周期的树行间杂草共同构成果园复杂环境。

1.1   图像灰度化

通过对果树不同时期分析可知，生长期果树多以绿色为主，其中掺杂红色果实，地面以黄色土壤与绿色杂草为主，其中掺杂掉落的红色果实；非生长期果树树干为灰褐色且顶部细小分支较多，地面以黄色土壤、落叶以及绿色杂草为主，不同时期果园颜色特征不同，但其共同点为背景天空颜色不变，故以红−绿−蓝(Red-green-blue，RGB)颜色空间蓝色分量(B分量)为特征进行天空与果树的分离。

在图像处理前，以图像左上角为坐标原点，水平向右为X轴正向，竖直向下为Y轴正向，沿Y轴将图像均匀分为四等份，由于图像两侧近视野区域不包含天空信息，属于无效处理区域，故去除左右各四分之一区域，只保留中间二分之一有效区域，灰度化效果如图1所示。

图  1  生长期与非生长期B分量灰度化效果

Figure  1.  Graying effect of B component in growth period and non-growth period

下载: 全尺寸图片 幻灯片

1.2   自适应阈值二值化

目前常用的二值化方法为最大类间方差法(Otsu算法)^[13-15]，其基本思想是根据图像灰度信息将图像分为目标与背景两类，以两者灰度方差达到最大时的灰度级作为分割阈值。由于非生长期时果树树干顶部细小分支较多，其灰度值与背景相近，若直接使用Otsu算法进行二值化处理会导致部分枝干信息丢失，出现错分的现象，不利于树行直线提取。针对这一现象，本文提出一种基于最大类间方差法的自适应阈值分割方法，其原理为：

假设存在阈值 $ {r_0} $ ， $ {r_0} $ 属于图像灰度 $ 0,1,2, \cdots, L - 1 $ 范围内(L为图像的总灰度级)，该阈值可将图像分为目标 $ {T_1} $ (小于 $ {r_0} $ 的像素点集合)与背景 $ {T_2} $ (大于 $ {r_0} $ 的像素点集合)，图像总像素数为 $ N $ ，则灰度级为 $ i $ 的像素点总数目为 $ {n_i} $ ，其出现的概率 $ {p_i} $ 为：

像素被分为目标 $ {T_1} $ 与背景 $ {T_2} $ 的概率S分别为：

两类灰度均值分别为：

故整幅图像的灰度均值m为：

图像的类间方差表达式为：

将式(6)代入式(7)化简得：

令 $ i $ 在 $ [0,L - 1] $ 范围内遍历，当 $ {\sigma ^2} $ 取最大值时的灰度级 $ i $ 对应为Otsu算法最佳分割阈值 $ {r_0} $ 。取灰度范围 $ [{r_0},L - 1] $ ， $G = {p_i}(i = {r_0}, \cdots , L - 1)$ ，令 ${G_{\max }} = \max \{ {p_i}\} (i = {r_0}, \cdots ,L - 1)$ ，若 $ {p_{{r_0}}} = {G_{\max }} $ ，则自适应分割阈值 $ r = {r_0} $ ；若 $ {p_{{r_0}}} \ne {G_{\max }} $ ，则取 $ {G_{\max }} $ 所对应的灰度级为 $ {r_1} $ ，取灰度范围 $ [{r_0},{r_1}] $ ， ${G'} = {p_i}(i = {r_0}, \cdots ,{r_1})$ ，令 $G_{\min }' = \min \{ {p_i}\} (i = {r_0}, \cdots , {r_1})$ ，若 $ {p_{{r_0}}} = G_{\min }' $ ，则 $ r = {r_0} $ ；若 $ {p_{{r_1}}} = G_{\min }' $ ，则 $ r = {r_1} $ ；若 $ {p_{{r_0}}} \ne G_{\min }' $ 且 $ {p_{{r_1}}} \ne G_{\min }' $ ，则 $ r $ 为 $ G_{\min }' $ 所对应的灰度级，其具体步骤如图2所示：

图  2  自适应阈值选取示意图

Figure  2.  Schematic diagram of adaptive threshold selection

下载: 全尺寸图片 幻灯片

如图3所示，沿图3a中直线(示例行)绘制其像素列灰度曲线(图3b)，取该直线上一细小枝干像素点(示例点)，在灰度曲线中对应出该点灰度值，由图3b可知，若直接使用Otsu算法，该枝干像素点被错分为背景类，使用本文自适应分割阈值可将其正确分割为枝干。将图像中小于自适应分割阈值 $ r $ 的像素点设为255(白色)，大于等于r的像素点设为0(黑色)，二值化效果如图4所示。

图  3  单行像素及其B分量像素列灰度曲线

Figure  3.  Single row pixel and gray curve of its B-component pixel column

下载: 全尺寸图片 幻灯片

图  4  Otsu算法与自适应阈值法二值化效果

Figure  4.  Effect of binarization of Otsu algorithm and adaptive threshold method

下载: 全尺寸图片 幻灯片

2.   动态ROI区域选取与特征点提取

由于以天空为背景进行处理时噪声较少，对特征点影响较小，故无需进行降噪处理，此处使用3×3的卷积核对二值图像进行膨胀处理，以进一步强化枝干信息，其效果如图5a所示。为避免部分果树顶端缺失对特征点提取产生干扰、提高特征点提取精度，根据树行顶部分布规律，本文提出一种动态阈值“V形”感兴趣区域(Region of interest，ROI)和基于ROI区域法线方向的特征点提取方法。

图  5  膨胀效果及其垂直投影曲线

Figure  5.  Effect of expansion and its vertical projection curve

下载: 全尺寸图片 幻灯片

由图5可以看出，左右两侧树行顶部呈直线分布，根据垂直投影曲线(图5b)最小值可区分左右树行，且最小值两侧曲线整体趋势与树行顶部分布规律相似，故通过确定投影曲线左右两侧最大值为纵坐标，取两端点与最小值点相连绘制两条直线，以两条直线为基准，取其上下各60像素宽度区域作为ROI，其结果如图6所示。

图  6  生长期与非生长期ROI区域

Figure  6.  ROI area in growth period and non-growth period

下载: 全尺寸图片 幻灯片

由图6可知，ROI由“V形”上边界和下边界组成，以左侧树行为例，以“V形”上边界为基准，沿“V形”左翼即上边界左侧直线从最小值点开始向左上方逐个取点作为基准点，过基准点沿上边界左侧直线法线方向向下寻找特征点，将遍历到的第1个像素值为255的像素点作为特征点，同理寻找右侧树行特征点。树行尽头由于视野较远存在图像畸变，导致在后续处理过程中产生较大误差，故以最小值点为基准点的2条法线与下边界所围区域为特征点无效提取区域，此区域内不进行特征点的提取，特征点提取示意图如图7所示。

图  7  特征点寻找示意图

Figure  7.  Schematic diagram of feature point finding

下载: 全尺寸图片 幻灯片

3.   可行驶区域识别

稳健回归是一种可用于异常点的检测并在样本存在较多异常点时替代最小二乘法的算法，常见稳健回归模型有随机采样一致性RANSAC算法和泰尔–森稳健回归(Theil-Sen robustness regression)，RANSAC算法是一种非确定性算法，当异常点较多时其拟合效果较差，而泰尔–森稳健回归对异常点具有较强的稳健性，基于此，本文提出基于双重稳健回归的思想进行树顶直线的拟合，以提高可行驶区域识别精度，其基本思路为：首先利用泰尔–森稳健回归对异常点具有较强稳健性的特点对偏离程度较大的异常点进行剔除，降低RANSAC算法的随机性；然后使用RANSAC算法对泰尔–森稳健回归剔除后剩余样本点进行直线拟合，提高得到正确模型的概率。

3.1   泰尔–森稳健回归剔除干扰点

为了更好地拟合特征点，提高树顶直线拟合精度，本文引入泰尔−森稳健回归剔除干扰点。泰尔−森稳健回归是一个参数中值估计器，适用于泛化中值，尤其是对多维数据的异常点有很强的稳健性，其一般回归方程为：

式中， $ \alpha $ 为拟合直线的截距， $ \beta $ 为拟合直线的斜率， $ \varepsilon $ 为模型的随机误差。其中：

从特征点中每次选取2个样本点，计算其所确定的斜率，直至所有样本点对选取完毕，若特征点数为 $ n $ ，则样本点对、所得斜率和截距理论上个数均为 $ C_n^2 $ 个，取所有斜率和截距的中值作为泰尔−森稳健回归拟合直线斜率 $ {k_0} $ 和截距 $ {\alpha _0} $ ，得到初步顶部直线方程。

由于非生长期果树枝干细小杂乱，特征点分布规律性较差，对直线拟合产生一定干扰，导致拟合误差较大，针对干扰程度较大的特征点，提出一种基于泰尔−森稳健回归拟合直线的自适应阈值特征点剔除方法，其基本思路为：通过特征点的横坐标 $ {x_i} $ 得到与其所对应的泰尔−森稳健回归直线方程纵坐标值 $ {y_{iT}} $ ，计算其与特征点纵坐标 $ {y_i} $ 的差值并取绝对值 $ {y_{iD}} $ ，所有特征点差值绝对值的平均值为特征点剔除阈值T，若 $ {y_{iD}} > T $ ，则将该点剔除，可用公式表达为：

3.2   关键点寻找及树行直线拟合

使用RANSAC算法对剔除后的特征点进行拟合，得到树顶拟合直线的最终斜率 $ {k_1} $ 和截距 $ {\alpha _1} $ ，确定树顶直线方程，图8为干扰点剔除及树顶直线拟合效果。受生长期树冠过大、非生长期部分枝干信息丢失以及拍摄角度等因素影响，树底直线斜率 $ k $ 由基础值 $ - {k_1} $ 和补偿值两部分组成。关键点的横坐标 $ {x_0} $ 取单侧树顶直线两端点横坐标的平均值，纵坐标 $ {y_0} $ 由“V形”下边界纵坐标最大值 $ {y_{\max }} $ 和剔除阈值组成，具体公式如下：

图  8  全周期干扰点剔除及树顶直线拟合效果

Figure  8.  Interference point elimination and fitting effect of straight line at tree top in full cycle

下载: 全尺寸图片 幻灯片

其中 $ \theta $ 为摄像头与水平方向的夹角， $ {y}_{0生} $ 和 $ {y}_{0非} $ 为生长期和非生长期时关键点纵坐标， $ {\left| {{k_1}} \right|_{\max }} $ 和 $ {\left| {{k_1}} \right|_{\min }} $ 分别为左右两侧树顶直线斜率绝对值的最大、最小值， $ {T_{\max }} $ 和 $ {T_{\min }} $ 分别为特征点剔除阈值中的最大和最小值，通过关键点 $ ({x_0},{y_0}) $ 和斜率 $ k $ 可确定可行驶区域边缘直线方程，并使用最小二乘法进行拟合，以确定果树行间可行驶区域，拟合结果如图9所示。

图  9  生长期与非生长期可行驶区域识别效果

Figure  9.  Recognition results of drivable area in growth period and non-growth period

下载: 全尺寸图片 幻灯片

4.   试验与结果分析

试验设备采用AMD Ryzen 74800H，主频2.90 GHz，内存16 G，Windows 64位操作系统的计算机，可行驶区域识别程序在Anaconda(Spyder)集成开发环境下使用python 3.8编写完成，图10为图像采集设备，试验中所使用的农业机器人宽45 cm，长60 cm，相机离地1.3 m，为降低计算量、提高识别速度，统一使用分辨率为640像素×480像素的图像。

图  10  图像采集设备

Figure  10.  Image acquisition equipment

下载: 全尺寸图片 幻灯片

本文选取不同天气状况对基于Otsu算法改进的自适应阈值分割方法进行试验，图11为生长期与非生长期基于B分量的分割及识别效果，以图11中6张图片为例，对其Otsu算法与本文自适应分割阈值进行比较分析，其对比结果列于表1。通过图表分析可知，基于B分量进行目标与背景的分离并不受天气状况以及季节的影响，即使处于阴天状况下仍能实现较好的分割效果。改进后的分割阈值大于等于Otsu分割阈值，且分割阈值根据每张图像自身各灰度级出现概率而定，良好的自适应性使得枝叶分割更为饱满，尤其是对非生长期而言，使用本文阈值进行分割能够将与天空背景相近且细小的枝干精确分割，提高后续拟合的精度。通过图11区域识别效果可知，本文算法受天气、光照等因素影响较小，在不同时期、不同天气状况下仍能达到较好的识别效果。

图  11  不同天气状况下全周期分割及识别效果

Figure  11.  Effect of segmentation and recognition under different weather conditions in full cycle

下载: 全尺寸图片 幻灯片

本文对生长期和非生长期晴天、多云、阴天天气下的顺光、逆光工况进行试验，选取泰尔−森算法和RANSAC算法与本文算法识别效果进行对比。以人工识别可行驶区域边缘线为参考，当一幅图像左右任一边缘算法识别结果与人工识别结果偏差角大于5°时，认为算法识别无效。表2为3种算法的试验结果对比，由表2中数据可知，泰尔−森稳健回归算法的准确率高于RANSAC算法，这是由泰尔−森较RANSAC而言对异常点有较强稳健性的特点决定的，但当异常点较多时，仍会对泰尔−森算法的斜率与截距选取产生影响。由于本文引入双重稳健回归的概念，依赖泰尔−森算法的稳健性对异常点进行剔除，降低了后续RANSAC算法拟合的随机性，使得本文算法较前2种算法有更高的准确率和更好的稳定性。本文算法较泰尔−森算法和RANSAC算法平均偏差角度分别减小了8.28%和9.88%，标准差分别减少了6.25%和22.89%，准确率分别提高了4.64%和10.49%，虽平均运行时间较前2种算法有所增加，但仍能满足农业作业中实时性的要求。

表  1  不同天气状况下全周期分割阈值对比示例¹⁾

Table  1.  Comparison examples of segmentation thresholds under different weather conditions in full cycle

图像序号 Image serial number	Otsu分割阈值( $ {r_0} $) Otsu segmentation threshold	$ [{r_0},L - 1] $	r1	$ [{r_0},{r_1}] $	本文自适应分割阈值(r) Adaptive segmentation threshold in this paper
1	143	[143, 255]	231	[143, 231]	217
2	132	[132, 255]	192	[132, 192]	173
3	131	[131, 255]	208	[131, 208]	141
4	143	[143, 255]	210	[143, 210]	158
5	158	[158, 255]	231	[158, 231]	171
6	115	[115, 255]	191	[115, 191]	115
1) ${r_1} = {{\rm{arc}}} {G_{\max } }$； $r = {\rm{arc}}{G_{\min } }$

下载: 导出CSV 

| 显示表格

表  2  3种算法的可行驶区域边缘拟合结果

Table  2.  Edge fitting results of drivable area for three algorithms

算法 Algorithm	图像帧数 Image frames	错误帧数 Error frames	平均偏差角度/(°) Average deviation angle	标准差/(°) Standard deviation	准确率/% Accuracy	平均耗时/ms Average time
泰尔−森算法 Theil-Sen algorithm	480	45	1.69	1.76	90.6	79.5
RANSAC算法 RANSAC algorithm	480	68	1.72	2.14	85.8	76.0
本文算法 Algorithm in this paper	480	25	1.55	1.65	94.8	82.4

下载: 导出CSV 

| 显示表格

通过不同的性能评价指标对本文算法的准确性进行详细分析，表3和表4分别统计了本文算法在每种工况下对随机120帧图像的准确率、平均耗时、平均偏差以及标准差。由表3和表4数据可知，不同时期逆光时的准确率均高于顺光时准确率，非生长期顺光时准确率最低且平均偏差角度和标准差最大，算法稳定性较差，逆光时准确率最高且角度标准差最小，虽平均偏差角度非4种工况中最小，但算法稳定性最好，除生长期顺光时平均偏差角度最小外，生长期顺、逆光时的准确率、平均偏差角度和标准差均处于中等水平。不同时期顺光时关键点平均偏差半径均小于逆光时平均偏差半径，生长期顺光时的关键点标准差最小，稳定性最好，逆光时标准差最大，稳定性较差，从整体来看，偏差角度与偏差半径误差均能满足农业机器人作业时精度要求。

表  3  不同工况下本文算法的拟合结果

Table  3.  The fitting results of this algorithm under different working conditions

时期 Period	工况 Working condition	图像帧数 Image frames	错误帧数 Error frames	准确率/% Accuracy	平均耗时/ms Average time
生长期 Growth period	顺光 Frontlighting	120	7	94.2	82.3
	逆光 Backlighting	120	5	95.8	82.4
非生长期 Non-growth period	顺光 Frontlighting	120	9	92.5	82.3
	逆光 Backlighting	120	4	96.7	82.5

下载: 导出CSV 

| 显示表格

表  4  不同工况下本文算法的性能评价指标

Table  4.  The performance evaluation index of this algorithm under different working conditions

时期 Period	工况 Working condition	图像帧数 Image frames	平均偏差角度/(°) Average deviation angle	标准差/(°) Standard deviation	关键点平均偏差半径/像素 Average deviation radius at key points	关键点标准差/像素 Standard deviation at key points
生长期 Growth period	顺光 Frontlighting	120	1.30	1.53	1.23	2.09
	逆光 Backlighting	120	1.35	1.68	2.69	3.63
非生长期 Non-growth period	顺光 Frontlighting	120	2.21	1.94	1.69	3.13
	逆光 Backlighting	120	1.35	1.15	2.53	2.65

下载: 导出CSV 

| 显示表格

5.   结论

本文基于以天空为背景的方法避免了树冠阴影、杂草以及地面其他干扰因素对识别精度和特征点提取效果的影响，利用背景天空特征B分量进行树顶与天空的分离不受天气状况及季节的影响，避免了算法对枝叶颜色的依赖性，使得该算法在阴天、生长期枝叶稀疏以及非生长期时同样适用，改进Otsu算法良好的自适应性使得被错分的细小枝干成功从背景中分离出来，提高了分割精度，解决了文献[8]中对树冠颜色及饱满度要求较高的问题。根据垂直投影设定“V形”ROI区域并进行特征点的提取有效地保留了树顶分布规律，提高了特征点提取准确率。

引入双重稳健回归的概念进行树顶直线的拟合，利用泰尔−森稳健回归对异常点较强的稳健性进行干扰点剔除，提高了后续RANSAC算法得到正确模型的概率和拟合精度，通过斜率变换关系得到可行驶区域边缘直线斜率信息，并进行关键点的寻找，使用最小二乘法拟合获得可行驶区域边缘直线。对不同生长期、不同工况进行试验分析，结果表明本文算法准确率高达94.8%，平均处理时间仅为82.4 ms，满足农业机器人果园作业中准确性与实时性的要求。