基于特征融合的果园非结构化道路识别方法

张彦斐 封子晗 张嘉恒 宫金良 兰玉彬

(1.山东理工大学农业工程与食品科学学院, 淄博 255000; 2.山东理工大学机械工程学院, 淄博 255000)

0 引言

随着农业现代化的推进[1],我国农业快速发展,智慧农业加速了农业机器人代替传统农业机械的步伐,促进了农业机械生产过程中自动化、智能化水平的提升[2],在减轻人工劳动强度的同时,还能降低作业成本、提高作业效率与精度[3]。道路的精准识别作为农业机器人实现自主导航、路径规划以及实施精细化作业的基础,对于农业机械在非结构化道路环境下是否能够安全可靠的运行具有重要意义[4]。

非结构化道路由于没有清晰的道路边界且易受阴影和水迹等干扰因素影响[5],导致道路区域和非道路区域难以区分。目前常见的道路识别方法主要依靠视觉传感器、激光雷达及多传感器数据融合等,由于视觉传感器获取信息丰富,可从图像中提取颜色、形状和纹理等特征以识别道路,不受多传感器协同作业时的干扰,且具备环境适应性强、误差累积小、性价比高以及细节识别精度高等特点[6-10],已成为农业车辆自主导航系统的研究热点。

在基于机器视觉识别道路研究中,文献[11-12]通过提取多种道路特征并结合支持向量机、随机森林等实现道路识别,该方法对结构化道路识别效果较好,但当面对道路环境相对复杂的非结构化道路时,识别效果并不理想。文献[13-16]通过道路形状特征、各颜色空间特征信息等进行道路特征的提取,从而限定路面的不可通行区域达到可通行区域检测的目的,并结合各种曲线模型、拟合算法完成对道路的识别。以上方法虽对非结构化道路检测效果较好,但对图像颜色依赖性较大,当道路与背景差异不显著以及阴影覆盖面积较大时算法检测失效。

针对果园环境中非结构化道路没有明显道路边界且周围环境复杂,单一使用颜色特征无法将道路和与其特征相似的干扰因素较好分离等问题,本文提出一种将颜色特征与纹理特征相结合的非结构化道路边缘识别方法,将全局特征提取与局部特征提取相结合,保证边界提取的精度与完整度,并使用基于距离和位置双重约束的两级伪特征点剔除方法对干扰点进行剔除,最后使用分段三次样条插值函数对边缘点进行拟合得到道路边缘曲线。

1 图像预处理

1.1 畸变矫正

使用RMONCAM林柏视-G600型相机采集分辨率为640像素×480像素的图像,此时图像位于像素坐标系中且存在一定畸变,为了得到畸变相对较小的图像,需进行相机标定确定相机内部参数以得到像素坐标系与世界坐标系的映射关系,并获取畸变参数以矫正由于镜头存在光学畸变造成的图像径向与切向畸变。

本文使用OpenCV中的相机标定方法[17]对20幅不同角度的11×8棋盘格图像进行标定,标定过程如图1所示,相机内部参数与畸变参数见表1,反向投影误差为0.025像素,对图像进行畸变矫正与变换缩放,效果如图2所示。

图1 相机标定过程

图2 图像畸变矫正

表1 相机标定结果

1.2 图像灰度化

矫正后图像位于RGB颜色空间中,对于非结构化道路而言,其周围环境较为复杂,除道路两侧植被外,还有阴影、土壤以及沙石等因素的干扰,而RGB颜色空间3个分量与亮度高度相关,易受光照影响,使用单一颜色特征难以将道路与非道路区域进行区分,故将RGB图像转换到更适合进行图像处理的HSV颜色空间[18]。

HSV颜色空间由色调(Hue,H)、饱和度(Saturation,S)和亮度(Value,V)3个分量构成,对图2b在HSV空间各分量的特征分布进行分析,图3为各分量的特征空间。

图3 HSV颜色空间各分量空间特征图

由图3可知,以图像左上角为原点,向右为x轴正方向,向下为y轴正方向,道路图像HSV颜色空间中H、V分量道路区域与周围植被区域灰度渐进,无明显灰度分级。灰度等级计算公式为

(1)

式中GH、GM、GL——HSV空间中各分量灰度分级后的高、中、低灰度级

Gmax、Gmin——HSV空间中各分量的最大、最小灰度

[]——取整符号

若将图像灰度等级分为高(GH)、中(GM)、低(GL)3级,H分量中道路与植被均处于GL灰度级(灰度范围为[0,59)),V分量中道路与植被均处于GH灰度级(灰度范围为[170,255]),不易区分道路与非道路区域。相较而言,S分量中道路区域主要位于GL灰度级(灰度范围为[0,85)),邻近区域主要位于GM与GH灰度级(灰度范围分别为[85,170)和[170,255]),两区域灰度分级较明显,利于后续道路边缘提取。对图4中不同工况和季节下图像中道路及其邻近区域进行灰度分级,其分级结果见表2。由表2可知,不同环境下S分量较H、V分量而言,其道路及其邻近区域均处于不同的灰度等级,且S分量受光照干扰程度较小,故本文选用HSV颜色空间中的S分量对图像进行灰度化处理,图2b灰度化效果如图5所示。

图4 不同环境下图像

图5 S分量灰度化效果

表2 不同环境下图像灰度分级结果

1.3 感兴趣区域选取

感兴趣区域(Region of interest,ROI)的选取要在满足减小计算量需求的同时,尽可能保留图像中的有效信息,以便后续图像处理。为能恰当选取ROI分割位置,本文提出一种基于滤波与梯度统计相结合的感兴趣区域动态提取方法。

由图3b可知,虽然道路与非道路区域灰度分级较明显,但非道路区域灰度变换剧烈且道路区域存在细小噪声。为方便后续梯度统计,使用双边滤波[19]对图像进行平滑处理,得到滤波后输出的像素值ID(i,j)为

(2)

式中 (i,j)、(u,v)——核内邻域像素与中心像素坐标

f(i,j)、f(u,v)——核内邻域像素与中心像素灰度

σs、σv——空间域核与值域核的标准差

D——以(u,v)为中心的滤波核窗口

双边滤波可保证在对图像消除低频分量噪声的同时保留边缘细节。滤波后,采用Sobel算子沿水平方向提取图像梯度,图6为滤波效果及图像梯度空间分布特征。

图6 双边滤波效果及Sobel水平梯度空间特征图

由图6b可知,因道路区域灰度变换均匀,该区域梯度基本为零,道路两侧及边缘存在梯度变换,且单行像素中道路区域所占比例由近视野区至远视野区逐渐变小,若对道路区域统计每行梯度非零的像素个数,道路尽头所处位置梯度非零个数应为最大。由于道路尽头位于图像中间区域,为减小计算量,沿图像水平方向选取中间1/3区域自上而下统计每行梯度非零的像素个数并存入数组N_zero[],计算数组最大值及其对应索引位置,若存在多处索引,则取最大索引并确定分割位置

(3)

式中L——ROI分割位置对应索引

h——输入图像高度

I——数组N_zero[]最大值对应索引位置

图7为S分量ROI。

图7 S分量ROI

2 图像融合与道路边缘提取

2.1 多特征图像融合

2.1.1纹理特征提取

Gabor滤波器[20-21]可对局部特征进行提取并对光照变化具有较强的鲁棒性,能够根据特定纹理方向提取有效特征并去除颜色特征无法识别的阴影部分。由于道路本身同样存在纹理特征,若设置的纹理提取方向过多易造成全图像纹理特征提取,导致道路与非道路区域无法区分。为保证在突显道路边缘特征的同时尽量不提取路面纹理,根据道路与边缘区域纹理角度信息,构建2个尺度为7、方向分别为0°和90°的Gabor滤波核并对S分量ROI进行滤波得到纹理图像g(x,y)。

(4)

其中

(5)

(6)

式中 (x,y)——像素坐标位置

θ、σ、γ——滤波核的方向、标准差与长宽比

b、λ、ψ——滤波核的带宽、正弦函数的波长与相位偏移

图8a、8b为滤波后所得纹理图像,由于θ=0°的原纹理图像对比度较小,为使纹理特征更便于观察,图8a为θ=0°原纹理图像做一次灰度正规化增强处理后效果,图8c为θ=0°原纹理图像与θ=90°纹理图像相加所得的多方向纹理特征提取图像。

图8 纹理特征提取

2.1.2图像融合及后处理

当路边土壤、沙石等与道路饱和度相近的干扰因素覆盖道路时,单一使用S分量无法将其与道路分离,易造成边界的误判或漏判。为保证边界提取的精度与完整度,本文提出一种基于图像颜色与多方向纹理特征相结合的图像融合方法,该方法主要包括以下两个关键步骤:

(1)纹理特征掩膜提取。采用Sobel算子计算纹理图像水平方向梯度,依据梯度图像中边缘位置较其他位置而言其梯度绝对值更大的特点,将图像梯度由最大值至最小值平均分为5份,提取梯度位于上2/5与下2/5区域的像素点,其余像素点值设为0,所提取到的特征掩膜如图9a所示。

图9 图像融合及二值化

(2)最大值法图像融合。将提取到的纹理特征掩膜各像素点值与S分量ROI最大类间方差法[22-24]分割阈值进行比较,根据比较结果对掩膜像素点值进行重置。

(7)

式中V′M(i,j)——重置后纹理特征掩膜像素点值

VM(i,j)——重置前纹理特征掩膜像素点值

VS(i,j)——S分量ROI的像素点值

TS——S分量ROI的最大类间方差法分割阈值

重置后的掩膜图像与S分量ROI逐像素点比较,取两者中较大像素点值生成融合图像,如图9b所示。

F(i,j)=max{V′M(i,j),VS(i,j)}

(8)

式中F(i,j)——融合后图像像素点值

融合后图像将全局特征提取与局部特征提取相结合,在颜色特征提取基础边界的前提下,使用纹理特征对精细化边界进行填充,能够有效提高边缘检测精度。

使用最大类间方差法将融合图像二值化,对二值图像采用3×3模板进行开运算,采用最大联通区域法先将其面积小于100像素的白色孔洞像素值设为0(黑色),再将其面积小于200像素的黑色孔洞像素值设为255(白色),以此达到去除噪声的目的,图9c、9d为图像二值化及其去噪后效果。

2.2 道路边缘提取

2.2.1特征点寻找及剔除

若选取二值图像中包含边缘信息的某一列像素,可知道路边缘点应位于像素值突变位置。将去噪后二值图像沿水平方向平均分为40个子区域,取每个子区域中间列作为特征点提取列,将每列自上而下检索到的像素值由0(黑色)至255(白色)的突变点作为候补特征点。

候补特征点中存在许多干扰点,如单列存在多个突变点、道路区域大面积噪声未去除导致伪特征点被提取等情况,为精准剔除伪特征点保留边缘点,本文提出一种基于距离与位置双重约束的两级伪特征点剔除方法。基于距离约束的一级伪特征点剔除具体步骤如下:

(1)查找x值最小候补特征点所在提取列(以下简称为“首点列”)上存在几个特征点,若只有一个特征点,则保留该点(以下简称为“保留点”);若存在多个特征点,保留y值最大的特征点,其余点剔除。

(2)顺序遍历首点列后的每一提取列,若该列不存在特征点,继续遍历下一列;若该列存在一个特征点,则保留该点;若存在多个特征点,计算该列所有候补特征点与上一个保留点之间的距离(以下简称“点距”),并计算最小点距及其对应在该列候补点中的索引位置,若自上而下最小点距点为该列最后一个特征点,则保留该点,其余点剔除,否则计算其后所有候补特征点点距与最小点距之间的差值,若点距差值均大于给定阈值,则保留最小点距点,其余点剔除;若存在点距差值不大于设定阈值的候补特征点,则保留符合要求的候补特征点中y值最大的特征点,其余点剔除。

道路环境较简单时,可实现一级伪特征点剔除后的保留点均为道路边缘点,但当道路环境较为复杂时,一级伪特征点剔除后的保留点除道路真实边缘点外还包含少量较当列其他候补特征点而言虽距离边缘较近但非边缘点的干扰点,针对这部分干扰点,使用基于位置约束的二级伪特征点剔除方法进行剔除,具体步骤如下:

(1)使用K-means[25-26]方法对一级伪特征点剔除后的保留点进行聚类,并计算左右聚类中心,根据2个聚类中心x值大小区分左右特征点。

(2)将左侧保留点x值不大于其聚类中心x值的点集设为p1,保留点x值大于其聚类中心x值的点集设为p2;将右侧保留点x值小于其聚类中心x值的点集设为p3,保留点x值不小于其聚类中心x值的点集设为p4。

(3)将点集p1、p4中y值最大的保留点y1max和y4max,以及点集p2、p3中y值最小的保留点y2min和y3min设为基准点,图10为K-means聚类结果及基准点。

图10 K-means聚类结果及基准点

(4)判断基准点y1max和y3min是否为其所在类别中x值最小的特征点,若是,则保留点集p1、p3中所有特征点;若不是,剔除其所在类别中x值小于基准点y1max、y3min相应x值的保留点。

(5)判断基准点y2min和y4max是否为其所在类别中x值最大的特征点,若是,则保留点集p2、p4中所有特征点;若不是,剔除其所在类别中x值大于基准点y2min、y4max相应x值的保留点。

(6)经过步骤(4)与步骤(5)剔除后所剩的特征点为最终提取到的道路边缘点,特征点提取与剔除结果如图11所示。

图11 特征点提取与剔除

2.2.2道路边缘拟合

由于非结构化道路边缘无固有形状,若使用线性拟合、边缘模型等方法进行拟合,其拟合效果与实际道路边缘误差较大,故本文使用分段三次样条插值法分别对道路左右边缘进行拟合,在保证分段点间连接平滑的同时避免了高次函数导致的龙格现象,使拟合结果更贴合实际边缘。

以左侧边缘点为例,设给定的节点(边缘点)坐标为(xi,yi)(i=0,1,…,n),相邻2个节点构成一个分段区间,n+1个节点共构成n个分段区间,每个分段区间的三次样条函数表达式为

Si(x)=ai(x-xi)3+bi(x-xi)2+ci(x-xi)+di

(x∈[xi,xi+1],i=0,1,…,n-1)

(9)

式中ai、bi、ci、di——各分段区间三次样条函数各项系数

则n个分段区间对应的三次样条函数S(x)表达式为

(10)

同时,三次样条插值函数S(x)需满足以下4个插值条件:

(1)函数S(x)穿过所有节点,即Si(xi)=yi(i=0,1,…,n-1)且Sn-1(xn)=yn。

(2)函数S(x)在所有节点(除第1个节点和最后1个节点)处零阶连续,即Si(xi+1)=Si+1(xi+1)(i=0,1,…,n-2)。

(3)函数S(x)在所有节点(除第1个节点和最后1个节点)处一阶连续,即S′i(xi+1)=S′i+1(xi+1)(i=0,1,…,n-2)。

(4)函数S(x)在所有节点(除第1个节点和最后1个节点)处二阶连续,即S″i(xi+1)=S″i+1(xi+1)(i=0,1,…,n-2)。

由式(9)可得各分段区间三次样条函数Si(x)的一阶导函数与二阶导函数分别为

S′i(x)=3ai(x-xi)2+2bi(x-xi)+ci

(11)

S″i=6ai(x-xi)+2bi

(12)

根据条件(1)与式(9)可得

(13)

根据条件(2)与式(9)可得

(14)

根据条件(3)与式(11)可得

(15)

根据条件(4)与式(12)可得

6aihi+2bi-2bi+1=0

(16)

式(14)～(16)中,i=0,1,…,n-2;hi=xi+1-xi。

本文使用非节点边界法给出插值函数的边界条件。第1段函数S0(x)和第2段函数S1(x)在第2个节点(x1,y1)处三阶导数连续,同时要求倒数第2段函数Sn-2(x)和最后一段函数Sn-1(x)在倒数第2个节点(xn-1,yn-1)处三阶导数连续,由此可得S‴0(x1)=S‴1(x1)和S‴n-2(xn-1)=S‴n-1(xn-1),即

(17)

联立式(13)～(17)即可求解分段三次样条插值函数S(x),得到左侧道路边缘三次样条插值曲线,同理,可得右侧道路边缘三次样条插值曲线, 图12 为道路边缘拟合结果。

图12 道路边缘拟合结果

3 试验及结果分析

试验在位于山东省淄博市沂源县山东理工大学与山东中以现代智慧农业有限公司共建的标准化苹果园中进行,试验设备采用Intel(R)Core(TM)i5-7200U CPU,主频2.50 GHz,内存12 GB,Windows 64位操作系统的计算机,道路边缘识别程序在Anaconda(Spyder)集成开发环境下使用Python 3.8编写完成,图13为图像采集设备,图像采集时相机光心离地1.3 m。

图13 图像采集设备

在晴天、阴天(存在车辆干扰)、顺光、逆光、冬季晴天和雨雪天气6种工况条件下对本文提出的基于图像颜色与多方向纹理特征相结合的图像融合方法进行试验效果分析。在其他处理步骤一致的前提下,设置S分量与多方向纹理特征图像为对照组,提取3种特征图像中的边缘点,检测算法识别边缘点与人工标识边缘点之间的纵向偏差绝对值dv,并计算偏差率

(18)

式中Nv——ROI图像纵向像素数

rd——偏差率

图14为多工况条件下3种特征图像边缘点提取结果,表3为多工况条件下3种特征图像边缘点平均纵向偏差绝对值与平均偏差率。

图14 多工况条件下3种特征图像边缘点提取结果

表3 多工况条件下3种特征图像边缘点提取误差计算结果

由图14可知,在晴天、阴天、顺光、逆光4种工况下,基于S分量与融合图像提取到的边缘点较完整且误差较小,基于纹理图像提取的边缘点缺失严重且误差较大。在冬季晴天和雨雪天气工况下,3种特征图像提取的边缘点均较完整,基于S分量与纹理图像提取的边缘点误差略大于基于融合图像提取的边缘点误差。表3数据表明,在多工况条件下,S分量平均纵向偏差为1.94～3.00像素,平均偏差率为0.72%～1.54%,纹理图像平均纵向偏差为2.21～93.38像素,平均偏差率为0.74%～43.23%,融合图像平均纵向偏差为0.92～1.80像素,平均偏差率为0.29%～0.71%。S分量、纹理图像和融合图像在6种工况下的平均纵向偏差均值分别为2.43、39.71、1.36像素,平均偏差率均值分别为0.99%、18.02%和0.54%,S分量的平均纵向偏差均值与平均偏差率均值均居中,纹理图像的平均纵向偏差均值与平均偏差率均值均最大,融合图像的平均纵向偏差均值与平均偏差率均值均最小。融合图像在6种工况下的平均纵向偏差均值较S分量与纹理图像的平均纵向偏差均值分别减小44.03%和96.58%,融合图像的平均偏差率均值较S分量与纹理图像的平均偏差率均值分别减小了45.45%和97.00%。由此可见,本文提出的基于图像颜色与多方向纹理特征相结合的融合图像较单一特征图像的边缘点提取误差与偏差率更小、提取精度更高,且在阴天工况下,存在车辆干扰时,基于S分量与纹理图像的提取方法无法对车辆所在区域进行识别,而基于融合图像的提取方法可有效提取车辆所在区域边缘点,实现障碍物的躲避。

图15为常见拟合方法最小二乘法、随机采样一致性法(Random sample consensus,RANSAC)与本文分段三次样条插值法对边缘点的拟合效果。由图15可知,本文采用的分段三次样条插值法为曲线拟合,较最小二乘法与RANSAC法的直线拟合更贴合非结构化道路实际边缘曲线。当阴天工况下,路面存在车辆时,最小二乘法与RANSAC法无法对车辆所在区域边缘进行拟合,不利于农业车辆行驶过程中障碍物的躲避,而本文使用的三次样条插值曲线可实现对车辆所在区域边缘的准确拟合,有助于提高农业车辆行驶过程中的安全性。对6种工况各选取10幅图像,获取各边缘点所在列的实际道路边缘与算法识别道路边缘之间的纵向偏差绝对值,使用平均偏差、平均偏差率、偏差率标准差作为评价指标,对3种拟合方法的拟合效果进行分析,具体数据见表4与图16。

图15 多工况下不同拟合方法道路边缘拟合效果

图16 多工况下不同拟合方法的偏差率标准差

表4 多工况下不同拟合方法道路边缘拟合效果评价指标

由表4可知,3种拟合方法中,RANSAC法拟合道路边缘在6种工况下的平均偏差与平均偏差率均最大,其平均偏差范围为2.05～4.37像素,平均偏差均值为3.16像素,平均偏差率范围为1.03%～1.38%,平均偏差率均值为1.21%,雨雪天气时平均偏差及平均偏差率均最大;最小二乘法拟合道路边缘在6种工况下的平均偏差与平均偏差率均居中,其平均偏差范围为1.66～3.25像素,平均偏差均值为2.64像素,平均偏差率范围为0.84%～1.16%,平均偏差率均值为1.02%,雨雪天气时平均偏差最大,顺光工况下平均偏差率最大;本文分段三次样条插值法拟合道路边缘在6种工况下的平均偏差与平均偏差率均最小,其平均偏差范围为0.53～0.90像素,平均偏差均值为0.66像素,较RANSAC法与最小二乘法分别减小79.11%与75.00%,平均偏差率范围为0.17%～0.37%,平均偏差率均值为0.26%,较RANSAC法与最小二乘法分别降低78.51%与74.51%,晴天工况下平均偏差与平均偏差率均最大。

由图16可知,在同种工况条件下,使用分段三次样条插值法拟合的道路边缘较最小二乘法与RANSAC法,其偏差率标准差最小。6种工况下,RANSAC法拟合道路边缘的偏差率平均标准差最大,为0.31%,最小二乘法拟合道路边缘的偏差率平均标准差居中,为0.23%,本文分段三次样条插值法拟合道路边缘的偏差率平均标准差最小,为0.09%,较RANSAC法与最小二乘法分别减小70.97%和60.87%,表明本文分段三次样条插值法具有更好的稳定性。

由图16可知,本文算法雨雪天气时偏差率波动幅度最大,其标准差为0.15%,冬季晴天时偏差率波动幅度最小,其标准差为0.01%,顺光、阴天、晴天、逆光工况下偏差率波动幅度居中,其标准差分别为0.12%、0.10%、0.09%、0.07%。由此可知,算法在冬季晴天时最稳定,在雨雪天气下的标准差虽为6种工况中的最大值,但仍能满足农业机械作业过程中的精度要求。6种工况条件下,算法平均单幅图像处理耗时分别为93.4、90.0、83.5、93.2、88.9、90.3 ms,单帧图像平均处理时间为89.9 ms,满足农业机器人作业过程中的实时性要求。

4 结论

(1)针对果园中非结构化道路无明显边界且道路边缘存在阴影、土壤和沙石干扰等问题,提出一种颜色与纹理特征相结合的非结构化道路识别方法。根据滤波与梯度统计相结合的动态提取方法对HSV颜色空间S分量进行ROI选取,并在此基础上使用多方向纹理滤波器提取纹理特征,综合利用最大值图像融合法、基于距离与位置双重约束的两级伪特征点剔除方法、分段三次样条插值法对道路边缘进行拟合以实现非结构化道路识别。

(2)试验结果表明,6种工况条件下,融合图像提取的边缘点平均纵向偏差均值为1.36像素,平均偏差率均值为0.54%,较单一特征图像而言,其提取误差与偏差率更小、提取精度更高。分段三次样条插值法的平均偏差均值为0.66像素,平均偏差率均值为0.26%,偏差率平均标准差为0.09%,较最小二乘法与RANSAC法而言,拟合精度更高且具有更好的稳定性。当路面存在车辆时,基于融合图像的特征点提取方法可有效提取车辆所在区域边缘点,三次样条插值法可实现对车辆所在区域边缘的准确拟合,有利于农业机器人作业过程中障碍物的躲避,可提高其行驶的安全性。6种工况下,算法单帧图像平均处理时间为89.9 ms,满足农业机器人作业过程中的实时性要求。