基于卫星图像的平原植保作业区域检测研究

邵长超,王思捷

（河南科技学院信息工程学院,河南 新乡 453003）

植保无人机作业主要受非植保作业时间、能量消耗、路程、植保作业的覆盖率和遗漏率等约束条件影响[1].目前植保无人机航线主要靠人工规划,在进行植保作业区域,尤其是遇到不规则作业区域时容易产生较大误差,致使错误估算施药量,增加环境污染和植保作业时间.因此在植保作业前需要对作业区域感知,获取作业区域的形状和面积等信息,为作业规划提供二维数学计算环境[2].

目前,已有国内外学者围绕相关问题进行了研究.针对作业区域边缘问题,Umaselvi 等[3]基于聚类的彩色图像无监督分割方法,对卫星图像中的植被和市区进行分类.张明杰等[4]利用大飞机平台获取的高分辨率航空影像对平坦地区较规整的田埂进行研究,以提取田埂界线信息.针对区域面积计算问题,杨静静等[5]采用经典数学定理鞋带公式（Shoelace Formula）及计算机图形学中向量积法计算区域面积.孙艺哲等[6]针对小块农田和不规则农田使用基于改进后的Alpha Shapes 算法农机作业面积测量方法,结果误差控制在3.5%以内.

可以看到,目前研究大都建立在高分辨率的遥感影像基础之上,民用获取尚有难度,使用分辨率有限的民用卫星图像,可有效地降低作业规划的前期准备难度与成本.本文基于分辨率有限的民用卫星图像,围绕作业区域边缘检测和面积计算问题进行研究,为作业规划和导航提供研究基础.

1 研究方法设计

总体设计是获取到卫星地图后对图像进行处理,获得作业区域的二值图像,最后运用像素比例计算实际作业面积,组织结构见图1.研究具体过程为：首先利用高德地图[7]获取卫星图像,对获取到的卫星图像进行预处理;然后使用边缘检测算子检测作业区域边界,获得边缘检测图像后叠加掩膜与运算得到带有绿色域的图像,之后将绿色色域转换为背景色;再后对处理图像进行形态学处理与端点检测修补,连接断点以获得完整连通区域,对连通区域进行轮廓绘制和填充;最后将二值图像进行像素比例计算,依据颜色占比与实际图像总面积换算可获得作业区域的实际面积.

图1 组织结构Fig.1 The structure of organization

为模拟真实植保无人机作业规划软件,运用高德地图开放平台JS API[8]开发了多边形区域面积简易测量工具,用以与图像处理得到的面积进行拟合度对比.区域面积测量工具采用Shoelace 算法（也称高斯面积公式）,公式以n×2 的矩阵形式表示多边形上顺序排列的顶点,行列式的计算又存在错位,形如所系的“鞋带”而得名[9].其缺点是无法直接应用在区域内有弧度的节点上,即使方法改造后也会丧失一部分精度,加之为手工固定测量点,测量精度无法保证[10].

1.1 边缘检测

本文使用Python 环境下Open C V 库对图像进行处理.为了得到最佳的结果,获取图像后使用了多种边缘检测算子进行对比检测.主要有一阶Roberts Cross 算子、Sobel 算子等,二阶Laplacian 算子、Canny 算子等.Sobel 算子在边缘定位方面有较好的表现,缺点是容易产生过多的像素边缘[11];Roberts 算子虽然在定位边缘的准确度上也相对较好,但在抗噪声方面表现稍有逊色.经过验证对比最终选择了具有较高的检测准确率以及信噪比,且图像边缘的连续性和完整性较好的Canny 算子[12].Canny 算子的检测过程为[13-14]：

过程一：将彩色图像转换为灰度图像;

过程二：对图像进行高斯滤波去除图像噪声;

过程三：使用一阶偏导的有限差分计算梯度的幅值和方向;

过程四：对梯度幅值进行非极大值抑制,寻找像素点局部最大值;

过程五：双阈值算法检测和连接边缘.

通过对Canny 算子分析,Canny 算子图像最终得到的为灰度图像,无法区分检测到的边缘为田埂或者绿色作物,容易影响到边缘的精度.针对该问题,本文在Canny 算子的基础上对其做出了适应性改进,主要思想是将原始图像和Canny 结果叠加掩膜后进行and 运算,得到带有绿色域的图像,之后将图像模式转换为HSV 颜色模型,并将绿色域转换为黑色背景色,最后得到无绿色域作物的图像,见图2.其具体的实现过程为：

图2 原始图像与Canny 图像Fig.2 Original image and Canny image

过程一：将原始图像进行灰度处理;

过程二：为了减少噪声对边缘检测结果的影响,采用自适应滤波核对图像去除噪声;

过程三：使用一阶偏导的有限差分计算梯度的幅值和方向;

过程四：对梯度幅值进行非极大值抑制,寻找像素点局部最大值;

过程五：双阈值检测,阈值设定为最小,保留最多细节;

过程六：Canny 灰度图像与原始图像利用掩膜（mask）进行二进制and 运算;

过程七：图像转为HSV 格式,将图像绿色域颜色转换为背景色;

1.2 断点修补

得到色域转换后的灰度图像,发现存在边缘像素断点问题.为获得图像封闭区域,需要对图像进行断点检测和断点缺陷修补.为减少算法开销,本文先使用形态学闭运算进行断点修补,闭运算过程为图像先膨胀后腐蚀,其作用是填充物体内细小空洞,连接邻近物体和平滑边界[15].经过验证,本文方案使用矩形膨胀方法,结构元2×2,迭代次数2 时可以将大部分断点封闭.算法闭运算公式

形态学处理后大部分的断点已经封闭,但也暴露出其缺点,如图3 所示.

图3 闭运算连通Fig.3 Closed operation

当两条轮廓像素距离相距特别近时,在进行闭运算操作时容易使两个轮廓黏连在一起,对轮廓的提取会产生一定影响.所以要使断点连接只在轮廓端点进行,需要对其进行端点检测,结果见图4.

图4 端点检测Fig.4 Endpoint detection

由图4 可知,以P1 为中心的八邻域,P1 周边数值的和为1 时为端点,和为0 时为孤点.得到端点后,断点修补的方法将非常明确,只需要判断两端点的距离设置一定的阈值,小于等于该阈值即可将断点用直线连接.

1.3 轮廓绘制填充

植保作业区域断点修补后,形成完整连通区域,再通过轮廓提取获得可参考的植保作业区域,而进行轮廓填充的主要目的为统计二值图像颜色像素比例.以图5 为例,主要过程为首先使用cv.findContours（）函数提取全部图像轮廓;再将提取到的全部轮廓使用cv.contourArea（）函数进行面积计算,通过计算只保留并绘制出的最大面积轮廓,也就是植保作业区域;最后将最大轮廓进行填充并统计二值图像颜色像素颜色个数与比例.

图5 轮廓提取与填充Fig.5 Contour Extraction and Contour Filling

1.4 作业区域面积

通过对面积计算函数分析,使用opencv 提供的面积计算函数contourArea（）时得到的值并非轮廓内亮点的像素值.根据官方文档提供的信息可知,contourArea（）使用格林公式计算面积,其复杂度比求非0像素数低,在像素分辨率有限且相对较低时,两种算法会出现不同结果：对于大或正方形的轮廓时,误差最小;对于较小轮廓或者椭圆形等不规则轮廓时,误差会较大.像素轮廓见图6.

图6 5×4 像素轮廓Fig.65×4 Image boundary

验证图6 中的5×4 像素轮廓,contourArea（）算法取轮廓像素中心点计算像素数,其面积ContourArea=12,计算轮廓内非0 像素数面积NonZero=20.结果可知,两种算法会产生不同的结果且误差较大.为减少两者误差,本文使用求非0 像素数获得作业区域面积,主要思想是二值图像黑白两色像素占比与实际总面积的比例换算,求得植保作业区域实际面积.具体算法过程如下

过程一 求出图像的实际长宽,公式为

式（2）、（3）中：L 表示实际长度,H 表示实际高度,l 表示宽像素数、h 表示高像素数,d 表示DPI,2.54表示英寸与厘米的比例关系,b 表示地图图像比例尺.

过程二 计算植保作业区域实际面积,公式为

式（4）中：S 表示作业区域面积,a 表示轮廓非0 像素数量,lh 表示图像总像素数,LH 表示土地面积.

2 验证与分析

获得到轮廓填充图像后最重要的一步是计算图像面积.因为卫星地图图像面积精确程度直接影响到最后的植保作业区域面积精度,所以首先需要求得卫星图像长宽与真实土地的比例关系.高德地图提供了一种地图缩放级别（zoom）,在高德开放平台的开发文档中,地图显示缩放级别范围默认为3～18,取值范围3～18[16].因开发文档中没有给出缩放等级与实际距离的比,不过通过验证计算可以得知,图像分辨率在96 DPI 环境下,每减小一个缩放等级距离像素比增大一倍,计算公式如下

式（5）中：b 表示比例尺,z 表示缩放等级.

根据计算公式可得比例尺与地图中的缩放级别的映射关系表1.

表1 缩放等级与映射关系Tab.1 Zoom level and mapping relationship

2.1 验证结果

地图验证采取高德地图缩放等级zoom=18（验证结果表明缩放等级越小,误差越大）下的卫星图像作为实验数据,地图采集地点在河南省新乡市红旗区某村农田,采集的模拟作业区域数为4 块.地图图像采集完成后,使用本文建构的方案对作业区域轮廓和面积计算方法加以验证,结果见图7 至图10.

图7 区域1 验证结果Fig.7 Verification results of region 1

图8 区域2 验证结果Fig.8 Verification results of region 2

图9 区域3 验证结果Fig.9 Verification results of region 3

图10 区域4 验证结果Fig.10 Verification results of region 4

以上验证的4 块区域图像中a 表示原始卫星图像,b 表示模拟作业区域轮廓提取,b 表示经过填充的区域,d 表示程序测量的结果.测量面积对比见表2.

表2 测量面积对比Tab.2 Results comparison

验证结果表明,本文方案对卫星图像处理后获得的轮廓信息,一定程度上抑制了边缘检测中的伪边缘,有效地清除掉了农田中的非种植区域,显示出农田真实轮廓.通过对表2 计算面积与程序测算面积的比对验证,两种方法在以上验证地块土地中的拟合程度最高达到了99.53%,从均值来看两种方式拟合度均值达到了95.78%.由本文方案得到的农田轮廓与面积信息表明具有较高的精确度,可以为植保作业提供数据参考.

2.2 误差分析

从验证结果上看,本文方案所得到的结果与手工测算具有些微误差,在此对误差产生原因进行分析：

（1）由于面积计算程序采用的算法产生的误差问题.由于鞋带公式的原理决定了无法避免测量误差,加之区域节点为手工添加也是图像拟合度降低的一个原因.

（2）卫星遥感影像微分纠正误差问题.卫星遥感影像偏斜是指由于地球自转,造成扫描行西移,使卫星影像呈平行四边形的偏斜现象[17].数字微分纠正指根据有关的参数与数字地面模型,利用相应的构像方程式,或按一定的数学模型用控制点解算,从原始非正射投影的数字影像获取正射影像[18].比较有趣的新闻就是“三峡歪了”[19],该新闻描述的内容就是微分纠正倾误差时出现了问题.本文方案利用高德获取到的地图都是经过微分纠正算法纠正过的图像.如果要进一步减少误差,可使用农田检测无人机提前获取作业区域的正视投影,但这无疑增加了部署时间,浪费人力物力.

（3）获取到的卫星图像分辨率低与比例尺误差的影响.由于获取到的卫星分辨率有限,图像内部分细节丢失,加上民用地图的比例尺精度要求有限等因素,所以无法避免误差.

（4）卫星地图投影产生的失真问题.高德地图使用的Web Mercator 投影,也是误差产生的重要原因.

3 小结

本文主要将提取的卫星地图进行图像处理,采用改进的边缘检测算子检测了植保作业区域边缘,对图像进行形态学处理与端点检测修补后得到封闭作业区域,最后基于像素比例计算得到作业区域的面积数据.经过对本文方法验证,使用分辨率有限的民用卫星图像,其结果也可较好地得到植保作业区域的轮廓范围,能一定程度的减少植保作业前期准备工作.再通过对本文方法与面积计算工具的对比,结果得到两者具有较高的拟合度,可以为植保无人机的作业规划提供较为精确的二维数学计算环境.