基于阈值分割的松材线虫病树计数与定位

2021-12-30余博文

武汉工程大学学报 2021年6期

余博文，黄巍

武汉工程大学计算机科学与工程学院，湖北武汉430205

松材线虫是我国危害较大的外来入侵物种之一，松材线虫病（pine wilt disease，PWD）是该物种造成的一种毁灭性森林虫害，通过松墨天牛（monochamusalternatus）等媒介昆虫传播于松树体内，进而引发松林病害。最新公告［1］显示目前松材线虫病已在我国18个省份（666个区、县、市）发现，导致大量松树枯黄死亡，造成严重经济和生态损失。不同密度的松材线虫病树处理方式并不相同，对于10%以上病树率的区域，需要组织专业队伍全面砍伐，对于病树密度不大、只是零星分散着病树的区域，可以只处理病树［2］。因此对区域病树的计数和定位是松材线虫病树检测与防治的关键。

随着无人机技术和计算机视觉技术的发展，利用无人机航拍图像和计算机视觉技术处理松材线虫病树计数与定位取得了较大的进步［3］。松材线虫入侵松树后，外部具有显著症状：先是针叶失水，褪绿，继而变褐，而后整株枯死，针叶全呈红黄色［4］。利用该特征对病树进行计数与定位，有助于病树的处理，防止疫情扩散。

传统计数方法大致分为无监督学习和有监督学习两大类。无监督学习算法主要利用外观相似性［5］和运动相似性［6］解决计数问题，但前者无法处理遮挡密度较高的计数问题，后者无法处理静止图像的计数问题。有监督学习分为基于检测的计数方法（counting by detection）和基于回归的计数方法（counting by regression）。基于检测的方法［7］一般通过检测图像中每个独立的目标来进行计数，这类方法首先生成候选框（bounding boxes），然后通过非极大值抑制（non-maximun suppression，NMS）筛选候选框［8］得到病树的区域，虽然有较好的效果，但筛选过程需要大量的计算。也有研究提出通过局部关系推测目标数量的方法［9-11］，虽然省略了传统基于检测的方法中非极大值抑制筛选候选框的计算量，但算法时间复杂度仍然较高。假设背景单一、目标均匀分布的前提下，蒙特卡洛方法［12］和形态学分析［13］可以达到较好的检测效果，但假设不成立的情况下，这些方法计数效果都不理想。基于深度学习的检测方法，包括两阶段（two-stage）的fast-RCNN［14］和一阶段（one-stage）的YOLO［15］检测方法，有较高的检测精度，但计算量大且需要大量数据支持［16］。基于检测的计数方法简单易行，但对于高密度、高遮挡的复杂环境不太适用［17］。基于回归的计数方法无需找到目标的具体位置，主要是寻找图像特征与目标数量的映射关系。对低密度计数图像，可以使用背景消除［18］分割出前景像素，然后使用最小二乘法［19］或神经网络［20］获得参数模型；对高密度计数图像，可以使用灰度共生矩阵［21］、切比雪夫矩等［22］描述纹理特征，再使用神经网络［21］或SVM［23］学习纹理特征与目标数量的映射关系。结合前景像素和纹理特征的回归计数方法［24-26］虽然避免了目标检测的困难，但同时也丢失了目标空间信息。密度函数［27］概念的提出弥补了基于回归的计数方法的不足，这种方法是通过图像特征与图像计数目标密度图的映射关系来对目标计数，有效保留了计数目标的空间信息，进一步优化了计数精确度［28］。

松材线虫病区区域较大，建筑、岩石、土壤等干扰因素较多，这要求计数定位算法克服复杂场景下的干扰，以提高病木的计数定位精度。本文对密度函数计数方法进行改进，利用色调（hue）、饱和度（saturation）、明度（value）（简称HSV）颜色空间的阈值分割、形态学处理、超像素对背景消除以突出前景像素目标，减少计数定位的干扰因素，从而提高松材线虫病树的计数与定位精度。

1 基于阈值分割的计数方法

1.1 密度计数

传统密度计数方法在医学领域中取得了显著的效果［27］，这些方法一般利用归一化高斯核函数（normalized Gaussian kernel，NGK）构建真实密度函数：

式中假设I1,I2,…,IN是N幅训练图像表示均值点在P处，协方差矩阵为σ2E的归一化高斯核函数，一般σ值取计数目标半径。

使用线性函数将每个像素点的特征向量映射到密度函数值，即带参数的密度函数为：

每幅图像Ii的每个像素点p都可以用一个K维向量表示∈R K，ω∈R K是密度函数模型线性变换的参数向量就是指定参数ω时密度函数的估计。

根据经验风险最小化原理（empirical risk minimization，ERM）求参数ω，同时加入正则化项，以防止过拟合：

式中D是距离函数，它是训练密度函数与真实密度函数差异的一种度量，ωΤω是复杂度惩罚项，避免求得ω值过大，正数λ为正则化参数，控制正则化强度。

这种密度函数模型适用于背景像素均匀且单一的医学细胞图像，但在松材线虫病树区域这种复杂环境的图像中无法达到理想的计数精度。这是因为松材线虫病树区域存在梯田、土壤、建筑等大量干扰因素，如图1箭头所指区域，这些干扰会对密度计数结果产生负面的影响。

图1 密度计数干扰因素示意图Fig.1 Schematic diagram of interference factors in density counting

如图2所示，为了减弱干扰因素的影响，提高目标计数的准确度，可以首先使用HSV颜色空间阈值分割方法预处理原始采集的图像，HSV颜色空间的阈值分割又分为形态学处理分割和超像素颜色均值分割。HSV颜色空间的阈值分割能够降低复杂图像中的干扰因素的影响，提高松材线虫病树的计数与定位精度。完整算法的结构如图2所示，输入图像首先经过HSV颜色空间的阈值分割处理，得到过滤图像，再使用密度计数对过滤图像进行处理，输出计数结果和对应的密度图像。

图2 基于阈值分割的计数方法流程图Fig.2 Flow chart of counting method based on threshold segmentation

1.2 HV通道过滤

本文使用HSV颜色空间的阈值分割直接去除背景像素，分割完成后，像素坐标点(x,y)处的像素值为：

式中f(x,y)表示H通道或V通道的值，Y表示阈值。当存在2个以上灰度阈值时，式（4）变为：

由于图像中除前景像素中需要计数的目标病树外，背景像素中还存在大量健康树木、土壤、梯田、建筑等干扰因素，这些干扰因素和目标病树的中心点的H和V通道值的散点图如图3所示。

图3中红、绿、黄、蓝、紫点分别代表病树、健康树、土壤、梯田、建筑。其中健康树木由于色调与目标病树明显不同，可以利用HSV颜色空间的H通道设定阈值去除，在数据集中随机选取50个病树中心提取分析H通道数值，绘制成直方图，如图4（a）所示。

图3 HV通道散点图Fig.3 Scatter plot for H-channel and V-channel

图4 病树中心像素的通道直方图：（a）H通道，（b）V通道Fig.4 Histogram of central pixel of infected trees：（a）H-channel，（b）V-channel

从图3中还发现土壤、梯田、建筑这些干扰因素虽然H通道值近似于目标病树的H通道值，但V通道数值大多数高于目标病树，可以利用HSV颜色空间的V通道设定阈值去除这些干扰因素，随机选取50个病树中心提取分析V通道数值，绘制成直方图，如图4（b）所示。

从图4中可以看出，病树H通道数值基本上都分布在9到35之间，V通道数值基本上都分布在29到46之间，可以设定9到35为H通道阈值，设定29到46为V通道阈值，对图1进行H通道过滤操作和V通道过滤操作，结果如图5所示。

图5 过滤通道：（a）过滤H通道，（b）过滤HV通道Fig.5 Filtered channel：（a）filtered H-channel，（b）filtered HV-channel

对比图5（a）和图1，可以发现健康树已被过滤，但还剩下大量土壤、梯田、建筑等干扰因素。经过HV通道过滤后，如图5（b）所示，土壤、梯田、建筑等干扰因素已被过滤，剩下的就是松材线虫病树初步估计区域。

1.3 形态学处理

原始采集得到的图像经过H和V通道过滤后可以再进行形态学处理以进一步提升密度估计的图像质量。由图5（b）发现经过H和V通道过滤的图像已经过滤了大部分土壤、梯田等干扰元素，但仍然存在许多接近点像素的细小干扰像素，且目标区域也有部分像素被错误过滤而形成“孔洞”。由于过滤图像是一幅二值化图像，因此利用形态学图像处理可以改善过滤图像的质量。

使用半径为2个像素的圆形结构元素C对H和V过滤图像B进行开操作，达到消除细小干扰像素的目的，执行开操作后的图像L由式（6）得到：

式（6）中：Θ表示腐蚀操作，⊕表示膨胀操作。使用C对L进行闭操作，达到修复目标区域内“孔洞”的目的，执行闭操作后的图像N由式（7）得到：

图5（b）所示H和V通道过滤图像经过形态学操作后的结果如图6所示。

从图6可以发现，大部分细小干扰像素已被过滤，目标区域部分过滤掉的像素造成的“孔洞”减少，图像更加圆滑，目标病树区域更加清晰。

图6 使用形态学处理的过滤HV通道Fig.6 Filtered HV-channel with morphological processing

1.4 超像素颜色均值图生成

超像素颜色均值分割即对原图像进行超像素处理再进行H和V通道过滤。分析图5（a）中H和V通道过滤图可以发现，这些多余的细小干扰像素的主要来源是由于梯田、土壤、建筑等干扰区域与病树的H通道和V通道数值接近，缺失的目标区域部分“孔洞”的来源主要由于某些病树像素的H通道或V通道数值与病树相差较大而被过滤。针对这些原因，可以对原图像做超像素处理，再执行分割过滤操作。

为了有效地利用相邻像素的纹理、颜色、亮度等特征将像素分组，本文使用简单线性迭代聚类算法（simple linear iterative clustering，SLIC）［29］对原始图像进行超像素分割。然后对每一个超像素进行均值化处理，得到超像素颜色均值图。图1的局部图像的超像素均值图如图7所示。

图7 超像素颜色均值Fig.7 Average color of superpixels

从图7可以发现利用超像素方法实现的像素值修正效果较好，将像素成功划分为具有相似值的病树、健康树、土壤等区域。对超像素颜色均值图进行HSV颜色空间的H和V通道过滤，结果与直接进行H和V通道过滤对比如图8所示。

从图8（b）可以发现，由于H通道和V通道数值接近病树而被保留下来的多余的干扰像素已被去除，由于H通道或V通道数值与病树相差较大而被过滤的缺失的目标区域部分已被补齐，图像更加圆滑，目标病树区域更加清晰。

图8 过滤HV通道对比：（a）不使用超像素处理的过滤HV通道，（b）使用超像素处理的过滤HV通道Fig.8 Comparison of filtered HV-channel：（a）filtered HV-channel without superpixelsprocessing，（b）filtered HV-channel with superpixelsprocessing

2 实验部分

2.1 数据集获取

实验数据为使用大疆无人机采集的湖北省某县松材线虫疫区图像，照相机型号为FC6310，其参数信息如表1所示。

表1 大疆无人机FC6310航拍相机及相片参数表Τab.1 Parameters of FC6310 aerial camera and photo of DJIUAV

数据集包括92张图像，无人机的飞行高度为350 m。为优化计算速度，对数据集图片作了重采样，得到的图片分辨率统一为972×728。按照7∶3的比例将数据集分为65张训练图像和27张测试图像，每张图像都有对应的ground-truth人工标注文件。

2.2 评价指标

为了对比经典的密度计数算法与本文算法性能差异，采用平均绝对误差（mean absolute error，MAE）、均方根误差（root mean squared error，RMSE）和绝对误差的方差（variance of absolute error，VAE）作为评价指标，平均绝对误差M、均方根误差R和绝对误差的方差V的定义如下：

式中N是测试图像数量，zi是第i张图像中实际病木数量，ei是第i张图像中算法估计的病木数量。通常M值和R值越小说明模型预测准确度越高，V值越小说明模型鲁棒性越高。

为使评价指标更加直观，使用计数准确率作为第4个评价指标，其计算公式如下：

2.3 对比实验

为验证本文算法的有效性，选择了Lempitsky的密度计数算法（density learning，DL）［27］、使用形态学处理分割的密度计数法（density learning with morphological processing，DLMP）、使用超像素颜色均值分割的密度计数法（density learning with super pixel，DLSP）进行对比实验。其中DLMP和DLSP均属于使用HSV颜色空间的阈值分割的密度计数法。

不同算法的评价指标如表2所示。

表2 不同算法计数误差对比Τab.2 Comparison of counting errors of different algorithms

从表2中可以看出，DLMP和DLSP算法的M、R和V3个评估参数都优于Lempitsky的DL密度计数算法［27］，相对于DL算法，DLMP和DLSP算法的MAE分别下降了10.9和12.5，DLMP和DLSP算法的RMSE分别下降了8.4和10.5，DLMP和DLSP算法的VAE分别下降了56.0和77.6，DLMP和DLSP算法的计数准确率分别提高了23.2%和26.6%。由此可见，使用HSV颜色空间的阈值分割可以有效提升密度函数模型的计数精度。

图9展示了不同算法在多建筑、多土壤、密林区域计数效果的密度图。为了便于比较分析，图9截取了完整图像的局部细节，图9中红色箭头表示目标病树漏检或误检。

从图9中可以发现，Lempitsky的传统密度计数算法DL对于松材线虫病区的病树存在大量误检和漏检，本文DLMP和DLSP算法相较于DL算法，对于多建筑、土壤、岩石等干扰项的计数区域计数效果良好，且密度图中标注的病树位置准确，只存在少量病树的误检和漏检，这证明了DLMP和DLSP有能力处理复杂背景下的计数问题。

图9 密度图对比：（a）航拍图像，（b）病树真实位置，（c）DL算法密度图，（d）DLMP算法密度图，（e）DLSP算法密度图Fig.9 Comparison of density maps：（a）aerial imagery，（b）ground truth，（c）density map generated by DL，（d）density map generated by DLMP，（e）density map generated by DLSP