1 引言

农业机械智能自主作业可提高农业生产效率、降低对劳动力的依赖,是农业发展的重要方向。从视觉图像中可靠、准确地分割出机械可行驶区域是农业机械在田间自主作业的关键。目前视觉导航的环境感知算法有两类:1)基于双目视觉和激光雷达的环境感知算法^[1-2];2)基于单目视觉的环境感知算法。因为单目视觉摄像头搭建成本低于激光雷达,且相对于汽车,农业机械作业时运动速度较慢,所以可使用设备成本相对低廉的单目视觉算法来进行导航。可行驶区域根据田地中土壤、杂草和作物的纹理、颜色等特征来提取。如文献[ 3-5]中采用了“2G-R-B”算法与最大类间方差、自动阈值分割的组合算法实现了绿色作物行或垄行的提取。该方法在有作物覆盖的田地中具有较好的时效性和抗噪性,但对于没有作物覆盖的待播种农地,其土壤颜色相近,该方法的分割效果不理想。文献[ 6-7]中根据耕地机留在土地的划痕特征、杂草和作物纹理的尺度差,采用小波分析与其他特征提取算法的组合能有效地识别出导航路径。文献[ 8]提出了通过训练反向传播(BP)神经网络,从不同的杂草环境下提取绿色作物垄行的方法。该方法抗干扰能力强,能够正确识别垄行中的作物并以此为依据提取导航直线。文献[ 9-10]中提出了使用卷积神经网络来识别并定位图像中待检测对象的方法,但由于垄行贯穿于整幅图像,因此该类方法的训练样本标注困难。未种植作物的垄行由采光、干湿度、土质、开垄方式等因素决定其颜色和纹理特征,在纹理和颜色特征上呈现多样性,故基于特定纹理或颜色设计的可行驶区域分割算法泛化性不强。由于垄行通常较长,在实地采集的图像上观察可发现,同一垄行的尺度由近及远逐渐变小,所以可使用尺度空间^[11]的思想构建垄行尺度空间的特征集,并在特征集不同尺度的图片中提取对应的垄行特征,从而达到分割垄和沟的目的。

2 尺度空间表示

一条垄行由凸起的垄以及左右相邻的两条沟组成。调整视觉采集摄像头角度,使其聚焦于工作区的垄行,让机械即将驶过的垄行贯穿于整幅视觉图像,有助于正确地检测出当前作业区垄和沟的位置。由此可采集得到图1所示的视觉图像。观察图像可知,同一条垄行的特征包含于不同的尺度中 ,其尺度由近及远逐渐变小。垄行靠近摄像头的部分尺度最大,距摄像头远端的部分尺度最小,所以可构建一个尺度空间,在不同的尺度下提取该垄行的特征。目前图像多尺度特征提取有很多的方法与运用,如文献[ 12-13]中所采用的方法。本文通过构建高斯差分金字塔的方式来提取不同尺度的垄行特征。

图 1. 垄行图像

Fig. 1. Images of ridge row

下载图片 查看所有图片

2.1 尺度空间原理

多尺度分解的基本思想是将原始信号嵌入到一个单参数的派生信号族中。Lindeberg^[11]提出了尺度空间理论,并指出高斯核函数的线性和位移不变性决定了高斯核函数及其导数是唯一可以作为多尺度分解的核函数。使用其构造的多尺度空间中,粗尺度结构能以良好的方式与细尺度相关联。二维高斯核表达式为

G(x,y,σ)=12πσ2exp-x2+y22σ2,(1)

式中:σ为尺度参数;σ²为尺度级别,尺度级别随σ的增大而增大。设输入图像为I(x,y),它的尺度空间是由一组不同尺度参数的G(x,y,σ)分别与I(x,y)进行卷积运算后得到。若以L(x,y,σ)表示卷积后的图像,“*”表示卷积运算,则表达式如下

L(x,y,σ)=G(x,y,σ)*I(x,y)。(2)

当σ=0时,高斯核函数变为一个二维脉冲信号, 它与I(x,y)进行卷积运算得到原图像本身。当σ增大时,I(x,y)与二维高斯核的卷积运算使得图像的平滑度增高,以至于图像中尺度级别小于σ²代表细节部分的高频分量丢失,留下尺度级别大于σ²代表语义部分的低频分量。

若用一个常数乘法因子k分隔两个相邻尺度,则高斯核函数的导数可近似表示为

∂G∂σ≈G(x,y,kσ)-G(x,y,σ)kσ-σ。(3)

由于kσ-σ是常数,省略后并不影响尺度特征的表示,所以可以使用G(x,y,kσ)-G(x,y,σ)近似表示高斯核函数的导数,

G(x,y,kσ)-G(x,y,σ)≈(kσ-σ)×∂G∂σ≈∂G∂σ。(4)

于是Lowe^[14]提出可通过一个常数乘法因子k分隔的两个相邻尺度的差分计算近似得出仅包含σ²~(kσ)²尺度级别特征的图像,

D(x,y,σ)=[G(x,y,kσ)-G(x,y,σ)]*I(x,y)=L(x,y,kσ)-L(x,y,σ)。(5)

2.2 尺度空间创建

高斯差分金字塔^[14]是一种用于构建图像尺度空间的方法。垄行尺度空间的创建应用了高斯金字塔的结构,但构建金字塔的结构参数则根据后续特征提取算法的需求进行了调整。细分的尺度空间表示有利于尺度上具有连续性的垄行特征的提取,但尺度空间表示层数的增多会增加信息的冗余量,也会影响到计算的时效性,从而对计算机性能也提出了更高的要求。

垄行高斯差分金字塔的创建过程如下。

1)创建垄行高斯金字塔。该算法中,为了详细反映出垄行特征在不同尺度上的分布规律,金字塔的阶数和层数都比较多。如图2所示,垄行的高斯金字塔结构共包含7阶(Octave),其中每阶有8层(Interval)尺度空间表达。第一阶的第一层图像为原图像的灰度图I(x,y),尺度级别为0,且图片的像素需大于600×600。往后六阶的第一层图片,则由相邻上一阶的倒数第二层图像进行下采样得到。每一阶中,除第一层外的其他7层图像,则由每阶的第一层图像与一组不同尺度参数的高斯核卷积得到。该组高斯核卷积的尺度参数可表示为

σ(s)=σo(s-1),s=2,3,4,…,8,(6)

图 2. 高斯金字塔结构

Fig. 2. Structure of Gaussian pyramid

下载图片 查看所有图片

式中:o表示阶序数,o∈[1,2,3,…,7];s表示层序数,s∈[2,3,…,8]。若L_os(x,y)表示垄行高斯金字塔,则其第o阶中除第一层外的尺度空间表示为

Los(x,y)=Lo1*G[x,y,σ(s)],s=2,3,4,…,8。(7)

2)创建垄行高斯差分金字塔(DOG)。如图3所示,高斯金子塔中同阶相邻的两幅图片相减可得到仅含有该尺度级别特征的图像。该垄行的高斯金字塔共7阶,每阶有8层尺度空间表示,对应地可以得到共7阶每阶含有7层图像的高斯差分金字塔。垄行高斯差分金字塔中的特征图像D_o's'(x,y)可表示为

Do's'(x,y)=Los(x,y)-Lo(s+1)(x,y),(8)

式中:o'表示高斯差分金字塔的阶序数,o'∈[1,2,3,…,7];s'表示高斯差分金字塔的层序数,s'∈[1,2,3,…,7]。

图 3. 高斯差分算法

Fig. 3. Gaussian difference algorithm

下载图片 查看所有图片

高斯差分金字塔特征图像的数值比较小,只有将该图片归一化才能得到肉眼可识别的图像。将高斯金字塔中所有图片归一化并使用最近邻插值法将其还原为原图像分辨率后,且同阶的7层特征图排为一列,可得到具有相同分辨率的垄行高斯差分金字塔特征图像组,如图4所示。

图 4. 归一化的垄行多尺度特征图像组

Fig. 4. Multi-scale feature image group of normalized ridge row

下载图片 查看所有图片

其中,为了保证分割算法能适应不同像素的图像并能得到噪声较小的分割图,高斯金字塔模型的选择遵循两个原则:1)该模型能适应低分辨率的图像;2)模型构建的尺度空间能涵盖较大的尺度范围。

基于第一个原则,该模型的阶数设置为7阶,共可进行6次下采样操作。因为50万像素的照片长宽为800×600,该照片执行6次步长为2的下采样后,像素降为19×10,再进行下采样意义不大。

基于第二个原则,每阶设置为8层,并使用(6)式设置高斯核卷积的尺度参数。使用(6)式的方法搭建得高斯金字塔,第7阶第一层高斯图像的参数σ初始值为σ₀×6⁶,而使用Lowe^[14]的搭建方法,第7阶第一层高斯图像的参数σ初始值为σ₀×2⁶,故使用(6)式搭建的高斯金字塔可获得更大的尺度范围。Lowe给出的算法中σ₀=1.6,但由于(6)式计算的参数σ(s)大于Lowe的方法,故(6)式的初始值被设置为Lowe方法初始值的0.5,即σ₀=0.8。使用(6)式确定的高斯差分金字塔中,各层特征图D_o's'(x,y)的K值可由下式确定。

K(s')=s'-1s'。(9)

该高斯差分金字塔模型的参数虽然不是最优值,但该模型构建的高斯差分金字塔特征组能满足特垄行征提取与分割的需求。

3 特征提取与分割

由于垄行受光照、干湿度、开沟纹理等因数的影响,其在灰度图中呈现出明暗相间的条纹。在本文开发环境中,像素的灰度值越大,该像素显示越白。若垄行灰度图的“垄”区相对“沟”区平均灰度值较大,则分割后“垄”在图像中显示为白色条纹;若垄行灰度图的“垄”区相对“沟”区平均灰度值较小,则分割后“垄”在图像中显示为黑色条纹。

本文提出了一种垄行特征分段提取的算法,算法的具体步骤如下。

1)对高斯差分金字塔进行中心化、灰度饱和化处理。即对垄行高斯差分金字塔中的每一张图像D_o's'(x,y)使用下式分别进行计算,计算后的得到的各图像以 Do's'1(x,y)表示。

Do's'1(x,y)=a×Do's'(x,y)-E[Do's'(x,y)]max[Do's'(x,y)]-min[Do's'(x,y)],(10)

式中:a为饱和系数;E[D_o's'(x,y)]为D_o's'(x,y)的平均灰度值;max[D_o's'(x,y)]为D_o's'(x,y)中的最大灰度值;min[D_o's'(x,y)]为D_o's'(x,y)中的最小灰度值。

2)对 Do's'1(x,y)使用最邻近插值法进行上采样,使其分辨率与原图像I(x,y)的分辨率相同。上采样后的图像用 Do's'2(x,y)表示。

3)将每一张特征图像 Do's'2(x,y)平均分为N段,其中N为正整数。为方便叙述,在本文中N=3,即将每一张特征图像 Do's'2(x,y)都平均分为上、中、下三段。对所有特征图像的上段进行加权求和运算,可得到上段垄行的特征矩阵。对所有特征图像的中段进行加权求和运算,可得到中段垄行的特征矩阵。在高斯差分金字塔特征组通过分段加权求和后得到的垄行特征矩阵中,数值大于1的区域为“垄区”,小于0的区域为“沟区”,数值在[0,1]区间的为不确定区域,以灰色斑点的形式出现在分割结果中。

在(10)式的计算中,对于本文使用的图像,中心化处理可以使代表沟的备选区域对应的特征值接近0或变为负值;去均值化后的图片乘以饱和系数a可使代表“垄”的备选区域对应的特征值超过1,代表“沟”的备选区域对应的特征值小于-1。通过多次实验对比,a取5时,能得到较好的效果。

算法步骤3)中,需要确定各图像求和时的权值。由图4可知,需检测的垄行的特征分布于第二阶至第七阶的各图像中。一幅图像中由上往下,位于图像的最上端的垄行距离摄像头较远,尺度较小,故特征主要集中在第二阶至第五阶的图像中;位于图像中间的一段垄行是作业机械即将经过的区域,是垄行特征识别的重点区域,其尺度相对于远端更大,故特征主要集中在第三阶至第六阶的图像中;位于图像底部的一段垄行距作业机械最近,尺度最大,图像中以垄为主,该段垄行是与上一垄行识别周期重叠检测的部分,在该周期的垄行识别中不做检测。

图 5. 分割的垄行图像。(a)上段垄行特征矩阵值的分布;(b)中段垄行特征矩阵值的分布;(c)上段垄行分割图;(d)中段垄行分割图

Fig. 5. Images of divided ridge row.(a) Distribution of feature matrix value of upper ridge row; (b) distribution of feature matrix value of middle ridge row; (c) segmentation of upper ridge row; (d)segmentation of middle ridge row

下载图片 查看所有图片

上段垄行的权值调整分两步。首先观察图4中上段垄行特征分布情况,根据其特征分布情况设置 Do's'2(x,y)各阶图像上段的权值;然后根据分割结果中噪点的数量及分布进行手工微调修正。经以上两步可得上段各阶的权值,如表1中w_h(o')所示。其中,w_h(1)=0表示 Do's'2(x,y)中,阶序数o'=1的共7层特征图片上段的权值;其他各阶特征图像上段的权值以此类推。将 Do's'2(x,y)中所有特征图像的上段乘以各自的权值并求和可得垄行上段特征矩阵,其值的分布如图5(a)所示。将垄行特征矩阵中小于0的值的设置为0,大于1的值设置为1,可得分割图像如图5(c)所示。

同理可得, Do's'2(x,y)各阶中段的权值如表1中w_m(o')所示。将 Do's'2(x,y)中所有特征图像的中段部份乘以各自的权值并求和,所得特征矩阵如图5(b)所示。根据得到的垄行特征矩阵分割出垄和沟的位置,分割图如图5(d)所示。

4 算法检验与分析

4.1 图像分割实验

验证用的垄行图像样本共40张,样本收集自17块不同土质的农地中,样本分辨率为700×620。不同的地块,其土壤纹理、开沟方式、垄沟走向、光照强度均不同。图像处理在计算机上进行,配置为:CPU型号i7-8850H,内存8G,采用Octave开发。

为了验证算法有效性,采用基于粒子群的最大类间差( PosOtsu)算法^[15]和空间核广义模糊C均值聚类(KGFCM)算法 [16进行实验对比。PosOtsu算法和KGFCM算法使用文献中的默认参数,分类数为2类。同时,为验证该分割算法分段参数N的取值对算法的影响,该验证实验中分别使用了N=3和N=6进行分割。为简便起见,用MsSegN3算法和MsSegN6算法分别记为本文算法N=3和N=6两个版本的简写。MsSegN3各段的权值在表1中已给出;MsSegN6各段的权值由表2给出。由于图像底部的两段垄行不做检测,表2仅给出了位于图像上端4段的权值,由图像上段往下,各段权值分别表示为w_h1(o')、w_h2(o')、w_m1(o')、w_m2(o')。部分分割结果如图6所示。

图 6. 垄行分割效果对比。(a)原图;(b) MsSegN3算法;(c) MsSegN6算法;(d) PosOtsu算法;(e) KGFCM算法

Fig. 6. Segmentation effect comparison of ridge row . (a) Original images; (b) MsSegN3 algorithm; (c) MsSegN6 algorithm; (d) PosOtsu algorithm; (e) KGFCM algorithm

下载图片 查看所有图片

垄行为全局性的特征,图像分割算法对土壤局部纹理特征的抑制能力决定了垄行分割的有效性。土壤纹理局部特征以小分辨率的“斑点”形态出现在分割结果中。“斑点”的数量越多则分割结果中连通域的数量越多。在各样本不同算法分割结果中,以黑色为背景的八连通域的数量级y如图7所示。该数量级y为标量,计算公式如下

y=ln(Q),(11)

式中:Q为分割结果中八连通域的数量。

图 7. 各算法分割结果中八连通域的数量级

Fig. 7. Number of eight connected domain in the segmentation results of each algorithm

下载图片 查看所有图片

结合图6和图7可以看出,PosOtsu算法得到的八连通域较其他三种算法多。30%样本的PosOtsu算法分割结果中八连通域的数量Q超过140,可以看出,该算法不能有效抑制土壤局部纹理特征。KGFCM算法分割得到的连通域虽然较少,但由于垄行上段尺度较小的部分不能被该算法有效分割,如图6(e)第一个分割结果所示,并因此导致了4个样本的分割结果中八连通域的数量Q=1,此时y=0,如图7所示。经统计发现,共计14个样本的垄行上段尺度较小部分不能被KGFCM算法有效分割。与以上两种算法相比,本文算法能有效抑制土壤的局部纹理所产生的噪声,同时能有效分割出不同尺度中的垄行。

4.2 算法主要影响因素

图6(a)中第一个样本图片,由于其垄的颜色较沟的颜色深暗,分割后“垄”区被显示为黑色的条纹,故该算法未能实现语义分割。

使用MsSegN3算法,单样本从加载到算法结束平均用时为0.439 s。使用MsSegN6算法,单样本从加载到算法结束平均用时为0.658 s。分段数量N的不同仅改变了特征组中各图片相加时的权值。在编程中,当输入图片大小固定时,如果将每一张特征图的权值均制作成与之同分辨率的模板,即以不共享权值的方式编程,则N取不同值不会改变算法的计算量,但这种编程方式会导致算法所需的运行内存增大;如果特征组中各特征图共享同段的权值,则编程时需要将各图像分为N段并使用循环对图像各段进行计算,计算量会随N的增加而增加,但该编程方式能节省运行内存的开销。取不同的N值,各段权值随机设置,算法使用共享权值的方式进行编程,通过批处理对40个样本进行测试,运行时间如图8所示。

图 8. N取不同值的算法运行时间

Fig. 8. Algorithm running time with different values of N

下载图片 查看所有图片

经实验发现,影响该算法分割精度的因素主要包括以下两点:

1)垄区自然光照的影响。若垄行是南北的走向,在有明显的阳光照射时,由于太阳不同时刻、光照的角度不同,导致背光面的“垄坡”在图像中显示较暗,该算法会将其误判为“沟”区。进而导致算法分割出的“垄”区的位置朝向阳面一侧发生偏移,偏移量为背阳一侧“垄坡”的宽度。

2)垄行局部纹理、颜色跃变的影响。垄行表面局部干湿度跃变和局部凸凹均会导致对应图像的纹理、颜色在局部出现跃变,并在分割结果中以黑色或白色的噪点显示。该类噪点可通过调整特征提取算法中的参数N来改善。N的取值越大,即同一条垄沟分的段数越多,则每一段垄沟对应的尺度区间范围越小,在其高斯差分特征集中,每一段垄行的特征能更精确地提取。由图6对比可知,同一垄行,当参数N=6时,分割图像中的噪点数明显少于N=3时的分割图像中的噪点数。但N值越大,需要设置的权值越多。

5 结论

本文提出一种垄行多尺度分割算法。该算法首先利用高斯差分金字塔结构构建垄行多尺度特征集,再对多尺度特征集中的图像进行灰度值饱和化、中心化、上采样处理,然后根据不同尺度下垄行特征的分布情况,将处理后的特征集中每一张特征图像分段,并分别进行加权求和运算得到垄行特征矩阵,最后根据垄行特征矩阵的值的将图像分割为“垄区”和“沟区”。 与PosOtsu算法和KGFCM算法相比,本文算法能在有效抑制土壤局部纹理产生的噪声的同时分割出图像中各尺度的垄行。

本文算法能够将尺度上由近及远逐渐变小的垄行分割为垄和沟,是一种有效的分割算法。后期希望对垄行多尺度特征集的阶数和层数进行优化;此外,目前该算法的饱和系数以及各阶垄行特征图像的权值需要手工调整确定,参数调整经验化,因此可通过改进该算法,从而根据参数N自动生成权值表并自动调整饱和系数;同时通过改进算法以减少分割图片中的噪声。