CN102542560A - 一种水稻移栽后密度自动检测的方法 - Google Patents

Abstract

本发明提供了一种水稻移栽后密度自动检测方法，该方法对所采集的实时前下视水稻田间图像进行分割，并利用自适应腐蚀和膨胀操作计算图像中水稻穴数，进而计算得到该块稻田中的水稻的种植密度。该方法以表征水稻穴的连通域属性作为判断依据，实时地对水稻生长图像进行计算，检测结果准确率高，对后期相关的农事活动具有重要的指导意义。

Description

一种水稻移栽后密度自动检测的方法

技术领域

本发明属于数字图像处理和农业气象观测交叉领域，具体涉及到一种水稻移栽后密度的自动检测方法，即以田间拍摄的水稻图像序列为对象，从图像特征上检测水稻种植密度的方法。

背景技术

水稻是我国主要的粮食作物之一，在我国南方广泛种植。一直以来，对于水稻各个发育期的观测主要是通过人工观测，受观测员主观因素的影响较大；同时由于水稻种植地域广、生长周期长，利用人工观测显然不够经济。因此，通过每天所拍摄的稻田图像，借助图像处理的手段，对其生长发育进行检测显得十分必要。水稻种植密度是田间水稻检测的一个重要环节，准确识别水稻种植密度，可以便与后期相关的田间作业，它是农业气象观测的一个重要内容。

2009年刘海娟在硕士论文“基于图像处理的水稻成熟期密度检测”中提出了一种利用图像检测水稻冠层密度分布的方法，利用检测到的冠层密度实时对联合收割机工作过程的调整喂入量，以便减少水稻在收割过程中的损失；2008年高建峰发表论文“辅以ETM的MODIS水稻种植面积监视研究”探索了ETM辅助下新一代卫星传感器MODIS在更小空间尺度下监视水稻种植面积。2005年田翠玲在硕士论文“水稻长势遥感监测与产量预报复合建模研究”中利用NOAA卫星的可见光通道和近红外通道，构造植被指数(NDVI)，结合地面样区实测叶面积指数，进行水稻的叶面积指数的遥感监测与产量的预报。以上方法中田翠玲和高建峰是利用遥感数据对水稻的种植面积进行确定，利用时间序列中水稻的叶面积指数，但由于遥感图像分辨率低，且容易受到云层、云阴影和气溶胶等影响，获得的水稻种植面积误差较大，无法准确获得水稻的种植密度。而刘海娟在硕士论文中只对水稻成熟期水稻的稻穗密度进行估测，此方法仍不能用来进行确定水稻的种植密度。然而，在农业气象观测领域，主要通过利用稻田的水稻图像实时的对水稻种植密度进行较准确的计算，以便及时指导后期间苗、估产等农事活动，而且稻田中水稻受浮萍等影响背景复杂，光照变化剧烈，上述利用遥感图像等方法都不可行。

发明内容

本发明目的在于提供一种水稻种植密度的自动检测方法，该方法利用田间水稻数字图像准确地检测出水稻的种植密度。

一种水稻移栽后密度自动检测方法，具体为：

对水稻图像进行二值分割；

初始化R个半径不同的圆形状结构元素StrEroded；

利用R个圆形状结构元素StrEroded对分割得到的二值图像进行腐蚀得到R个腐蚀结果bwerode_tR；

创建L个边长不同的方形形状结构元素StrDilate；

利用L个方形形状结构元素StrDilate对腐蚀结果bwerode_tR进行膨胀，得到R×L个膨胀结果bwdilate_tRL；

检测各膨胀结果bwdilate_tRL的异常连通域数量；

选取异常连通域数量从急剧下降到缓慢下降的转折点对应的方形形状结构元素作为最优结构膨胀元素；

利用最优膨胀元素对腐蚀结果bwerode_tR进行膨胀；

从利用最优膨胀结构元素对腐蚀结果bwerode_tR进行膨胀的结果中选取连通域数量最多对应的圆形状结构元素为最优腐蚀结构元素；

利用最优腐蚀结构元素和最优膨胀结构元素对二值图像进行先腐蚀再膨胀得到BWED_t；

统计BWED_t的连通域个数RegionNum_t；

计算水稻密度 ${\mathrm{RiceNum}}_t={\mathrm{RegionNum}}_t\×b\×\frac{1}{0.0015\×\mathrm{SceneArea}},$ b为每个稻穴中植入水稻数量，SceneArea为相机场景面积。

进一步地，所述异常连通域是指满足公式(1)或(2)的连通域；

Distance≤K₁×MinDist    (1)

Distance为本连通域与其它连通域质心的最短欧式距离，MinDist为所有连通域的Distance的均值，0.3≤K₁≤0.6；

$\left\{\begin{array}{c}\mathrm{Pixnum}\≤K_2\×\mathrm{MinPixnum}\\ \mathrm{Pixnum}\≥K_3\×\mathrm{MinPixnum}\end{array}\right.---\left(2\right)$

Pixnum为本连通域的像素数，MinPixnum为所有连通域的Pixnum的均值，0.15≤K₂≤1，2≤K₃≤3。

进一步地，在计算水稻密度步骤之前对所述BWED_t作去噪或/和连通域断开操作。

本发明的技术效果体现在：本发明自动对所采集的实时前下视水稻田间图像进行分割，并利用自适应腐蚀和膨胀操作计算图像中水稻穴数，进而计算得到该块稻田中的水稻的种植密度。该方法以表征水稻穴的连通域属性作为判断依据，实时地对水稻生长图像进行计算，检测结果准确率高，对后期相关的农事活动具有重要的指导意义。

附图说明

图1是训练阶段的流程图；

图2是水稻叶面彩色图像分割的流程图；

图3是检测水稻种植密度的流程图；

图4是训练图像集中随机挑选的一张；

图5是图4中的图像人工阈值分割的结果；

图6是待检测序列中随机挑选的一张图像；

图7是图6中图像的最终分割结果；

图8是图6分割结果自适应膨胀腐蚀处理后的结果；

图9是图6分割结果腐蚀结构元素大小为1时膨胀检索曲线；

图10是相机观测场景示意图。

具体实施方式

下面结合附图来详细说明本发明的较佳实施例，本发明中使用图像序列采集与离地面高5米的相机，镜头焦距为12毫米，水平视场角为46度，垂直视场角为35度，相机分辨率不低于400万像素。实施例以每一天为一检测时段，每一检测时段内拍摄w张水稻图像(w＝5)。每天为一检测阶段，有利于识别水稻的主要关键生长期。此发明旨在检测水稻种植密度以及监视稻田的密度变化，此处的水稻的密度只需三天左右检测一次。整个方法分为训练阶段、分割阶段和检测阶段。

1.训练阶段：通过前一年拍摄的稻田水稻的历史图像数据，统计图像中水稻叶面的H(色度)颜色特征随着Y(亮度)的变化规律，其流程如图1所示，具体步骤如下：

(1)样本图像生成，从往年历史图像数据中，随机选取各种光照条件下拍摄的稻田水稻图像100幅左右，将原图的R-G-B颜色空间转换为r-g-b颜色空间，再将r与g两个通道转化到极坐标系内，在极坐标系中利用方向与半径两个参数进行阈值分割，将水稻叶面区域保留，其它区域置为纯白色，得到水稻叶面的样本图像。图4是训练图像集中随机挑选的一张，图5是此张图像的人工阈值分割的结果。

(2)生成水稻叶面Y-H颜色特征统计表，水稻叶面Y-H颜色特征统计表包含信息有：水稻叶面像素点在各亮度下对应的色调均值和色调方差。在生成水稻叶面样本图像之后，统计其中水稻叶面在H-S-V和Y-Cb-Cr颜色空间中的变化规律，其中Y表示亮度，Cb表示蓝色色度分量，Cr表示红色色度分量，H表示色调，S表示饱和度，V表示亮度，具体过程如下：

设样本图像集为S_m，m＝1…150，图像中的第i行与第j列像素点相对应的在Y-Cb-Cr颜色空间的色彩分量Y_m(i，j)、Cb_m(i，j)、Cr_m(i，j)和在H-S-V颜色空间的色彩分量H_m(i，j)、S_m(i，j)、V_m(i，j)，依次从每张样本图像的每个像素点查找，记num＝1，若H_m(i，j)≠0或S_m(i，j)≠0或V_m(i，j)≠1(背景区域为纯白色)，则获取此水稻叶面像素点特征为[Y_num，H_num]，且令num＝num+1。依次查找直到终止，将最后的num记为TotalNum。于是上述过程得到一个TotalNum×2的水稻叶面像素特征矩阵F^TotalNum×2。

为充分利用水稻叶面在各种光照条件下的颜色特征，将F^TotalNum×2矩阵按行进行随机排列，再截取其中的前RowNum行特征，记为f^RowNum×2。定义全零矩阵z^256×RowNum，依次从i＝0，2，…，255，在f^RowNum×2(也即[Y_RowNum，H_RowNum])中查找Y＝i的所有H_m，m＝1，2，…N_i，其中N_i表示f^RowNum×2中Y＝i的像素数目，最后令z^256×RowNum的第i行元素从左端开始按顺序定义为H_m，m＝1，2，…N_i，N_i＜RowNum，不足的位置取0。

在像素点亮度为Y＝i，i＝0，1，…，255的情况下，可以求出水稻叶面的色调分量H的均值和方差，公式分别为：

${\stackrel{\‾}{h}}_i=\left\{\begin{array}{cc}\frac{{\mathrm{\Σ}}_{j=1}^{\mathrm{RowNum}}{z}_{y,j}}{N_i},& N_i\≠0;\\ 0,& N_i=0;\end{array}\right.$

和

${\mathrm{\σ}}_i^2=\left\{\begin{array}{cc}\frac{{\mathrm{\Σ}}_{j=1}^{\mathrm{RowNum}}{(z_{y,j}-\stackrel{\‾}{h_y})}^2}{N_i},& N_i\≠0;\\ 0,& N_i=0;\end{array}\right.$

其中N_i表示RowNum中对应于Y＝i的像素点个数，RowNum的取值区间为[2000，5000]，此处我们选择RowNum＝2200，获得了较好的结果。

因此，可以生成水稻叶面的Y-H颜色特征统计表：

$\mathrm{THTable}=\{y,{\stackrel{\‾}{h}}_i,{\mathrm{\σ}}_i^2|i=\mathrm{0,1},\·\·\·,255\}.$

上述训练阶段完成后，我们已经得到水稻叶面颜色的先验信息记为：

2.分割阶段，利用训练阶段生成的水稻叶面Y-H颜色特征统计表，对待检测的实时前下视稻田水稻图像进行水稻叶面的分割，流程图如附图2所示，具体操作步骤如下：

(1)利用水稻叶面的Y-H颜色特征统计表，对实时的前下视棉田视频图像进行水稻叶面区域进行粗选。将待分割图像序列记为D_t，t＝1，2，…n，文中所有下标t均代表待测图像序列中的每张图像。图像中的第i行与第j列像素点相对应的在Y-Cb-Cr颜色空间的色彩分量Y_t(i，j)、Cb_t(i，j)、Cr_t(i，j)和在H-S-V颜色空间的色彩分量H_t(i，j)、S_t(i，j)、V_t(i，j)，定义行列数与S相同的全零矩阵逻辑矩阵SL_t。

令

deltaH_t(i，j)＝|H_t(i，j)-meanH_t(Y_t(i，j)，1)|若

deltaH_t(i，j)＞1.9×standardH_t(Y_t(i，j)，1)

则此像素点为水稻叶面背景像素点，取SL_t此位置记为1。由此得到粗选后水稻叶面区域的分割结果逻辑矩阵SL_t。

(2)将待分割图像D_t转化到L-a-b颜色空间，L表示亮度对应值，a表示红或绿对应值，b表示黄或蓝对应值。将每个像素点在L-a-b空间下的a值和b值作为此像素点的特征进行k-means聚类，指定为五类。获得原图S的聚类标签CLabel_t，CLabel_t是行列数与D_t相同的单通道矩阵，且其中每个元素为D_t中对应位置像素点k-means的聚类标签。

(3)将上一步CLabel_t中的每一类对应点的区域与粗选后的分割结果SL_t中的元素值为0的区域分别进行叠加，计算每一类区域与粗选结果的重合度MDegree_t(r)，r＝1，2，…5；定义全零矩阵逻辑矩阵BW_t，行列数与SL_t相同，且令

(4)连通域标记和去噪，利用标记算法，对二值图像BW_t进行8邻域连通域标记，第k个连通域记为BW_t(k)，连通域的面积即像素个数记为num_t(k)，根据如下的公式：

去除较小的连通域以降低浮萍等水生植物的干扰。面积阈值Threshold根据需要进行选择，一般为[8，20]。一般可取Threshold为10。操作完成后得到最后的BW_t。查找BW_t中(i，j)位置为0的点，取S中此位置像素点为黑色，于是我们得到图像D_t中水稻叶面的最后分割结果DResult_t，图像序列中随机选取的一张图像如图6，则它的分割结果如图7所示。

本实施例分割阶段是利用历史统计的Y-H颜色特征统计表进行水稻叶面区域提取，这是优选方式，但不局限该方式，现有的图像分割方法均可使用本发明。

3.检测阶段，由分割阶段我们获得二值图像BW_t，使水稻叶面对应像素点值为1，背景像素点值为0。将二值图像BW_t进行自适应形态学操作膨胀腐蚀，然后利用连通域标记自动计算相机场景中水稻的穴数。进而根据场景面积，以及每穴的水稻株数，计算每亩的水稻种植数量。

(1)创建圆形状结构元素StrEroded，使其半径为R，R＝1，2，…，5。利用StrEroded对BW_t进行腐蚀(van den Boomgard，R，and R…van Balen，″Methods for FastMorphological Image Transforms Using Bitmapped Images，″Computer Vision，Graphics，andImage Processing：Graphical Models and Image Processing，Vol.54，Number 3，pp.254-258，May1992.)，结果定义为bwerode_tR，R＝1，2，…，5。

(2)创建方形状结构元素StrDilate，使其边长为L，L＝1，2，…，20。利用StrDilate对bwerode_tR进行膨胀(van den Boomgard，R，and R.van Balen，″Methods for FastMorphological Image Transforms Using Bitmapped Images，″Computer Vision，Graphics，andImage Processing：Graphical Models and Image Processing，Vol.54，Number 3，pp.254-258，May1992.)，结果定义为bwdilate_tRL，R＝1，2，…，5，L＝1，2，…，20。利用上述过程每张水稻分割图像BW_t可以获取RL＝5×20张膨胀腐蚀结果，后续步骤将从中自动择优，即确定最适应腐蚀膨胀结构元素大小参数。

(3)对于每个bwdilate_tRL检索其连通域属性，设bwdilate_tRL中有Region个连通域。计算bwdilate_tRL中各连通域质心位置(x_i，y_i)，以及每个连通域的像素点数量Pixnum_i，i＝1，2，…，Region。计算各个连通域质心之间的欧氏距离，进而可以获取每个连通域质心与其它连通域质心的最短欧式距离Distance_i，i＝1，2，…，Region，以及Distance_i的均值MinDist。对每个i利用判决公式：

Distance_i≤K₁×MinDist，i＝1，2，…，Region，0.3≤K₁≤0.6(优选0.56)如UnusualDist个满足上述公式，则记bwdilate_tRL中质心距离异常连通域个数为UnusualDist。计算Pixnum_i，i＝1，2，…，Region的均值为MinPixnum。对每个i利用判决公式：

$\left\{\begin{array}{c}{\mathrm{Pixnum}}_i\≤K_2\×\mathrm{MinPixnum}\\ {\mathrm{Pixnum}}_i\≥K_3\×\mathrm{MinPixnum}\end{array}\right.,i=\mathrm{1,2},\·\·\·,\mathrm{Region},$

0.15≤K₂≤1(优选0.2)，2≤K₃≤3(优选2.5)

如UnusualPixnum个满足上述公式，则记bwdilate_tRL中质心距离异常连通域个数为UnusualPixnum。利用连通域质心距离特征与像素点数目特征，定义bwdilate_tRL中异常连通域数量为UnusualNum＝UnusualPixnum+UnusualDist，计算UnusualNum的差商并记为diff。腐蚀结构元素大小为1时不同膨胀结构元素大小下异常连通域曲线如图9所示，图9中最左端点显示异常连通域数目超过连通域总数，这是由于质心距离过近的异常连通域同样可能是面积过小的异常连通域，即文中上面对异常连通域总数的定义有关，上述定义的异常连通域总数有利于后面水稻密度的检测。

(4)由上叙述，对每个bwerode_tR，R＝1，2，…，5分别进行了20次膨胀，下面获取每个bwerode_tR的最优膨胀结构元素大小。最优膨胀结构元素要保证将腐蚀段的连通域合并在一起，而且保证不能过膨胀。随着结构元素的增大，检测到的连通域异常数量变化会开始的时候急剧下降，而后缓慢下降(如图9所示)，急剧下降部分为合并大量腐蚀断的连通域造成的，缓慢下降部分为过膨胀时造成，于是两个阶段的转折点对应于最优膨胀结构元素的大小。利用上部分中计算得到的数组diff，由：diff＞-15，判断diff数组中第一个满足上式的元素位置，于是将其次序的大小定义为bwerode_tR的最优膨胀元素大小strel_dilate_opt。对上述的bwerode_tR用大小为strel_dilate_opt的方形膨胀元素进行膨胀。然后重复利用步骤(3)的方法计算bwerode_tR用最优膨胀元素膨胀后的连通域总数Num_dilate_tR以及异常连通域数量为UnusualNum_dilate_tR，即对每张分割结果图BW_t，由上得到正常连通域数量为UsualNum_dilate_tR＝Num_dilate_tR-UnusualNum_dilate_tR，R＝1，2，…，5。

(5)对每张分割结果图BW_t下面获取其最优腐蚀结构元素大小。最优腐蚀结构元素要求保证可以不把BW_t中小连通域腐蚀掉，即保证接合最优的膨胀元素膨胀后可以得到尽可能多的连通域。于是利用(4)中获得的UsualNum_dilate_tR，显然UsualNum_dilate_tR，R＝1，2，…，5为正整数数组，判断UsualNum_dilate_tR数组中最大元素的元素位置，于是将其次序的大小定义为BW_t的最优腐蚀元素大小strel_erode_opt。

(6)利用(5)中获得的最优腐蚀元素大小strel_erode_opt，以及(4)中获得的最优膨胀元素大小strel_delate_opt。对分割结果图BW_t进行先腐蚀再膨胀得到BWED_t，此时处理结果图BWED_t中包含有少部分很小的连通域(噪声)需要去除，检索BWED_t中连通域像素点数量BWPixnum_i，i＝1，2，…，Region′，和BWPixnum_i的中值BWPixnum_med，其中Region′是连通域个数。将满足条件：

${\mathrm{BWPixnum}}_i<\frac{1}{3}\&times;\mathrm{BWPixnum}\_\mathrm{med},i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,{\mathrm{Region}}^{\&prime;}$

的连通域像素点对应值置为0(去掉此类连通域)，获得每张分割结果图BW_t的粗选自适应腐蚀膨胀结果BWEDR_t。

(7)粗选自适应腐蚀膨胀结果BWEDR_t中会少部分BW_t中距离相近的水稻穴位连通域由于膨胀而连接在一起，此处应用二值图像的欧式距离变换尽可能将连接在一起的连通域断开。二值图像的距离变换用来计算二值图像中每个值为0的像素点与值为1的最近邻像素点的距离(Breu，Heinz，Joseph Gil，David Kirkpatrick，andMichael Werman，″Linear Time Euclidean Distance Transform Algorithms，″IEEE Transactions onPattern Analysis and Machine Intelligence，Vol.17，No.5，May 1995，pp.529-533.)。将BWEDR_t值为0的像素点置为1，值为1的像素点置为0，然后应用二值图像的距离变换得到BWEDRS_t。计算BWEDRS_t连通域中像素点的距离ReigionDist_i，i＝1，2，…，Region″，以及距离的中值ReigionDist_med，其中Region″是连通域个数。将满足条件：

ReigionDist_i＜1.9×ReigionDist_med，i＝1，2，…，Region″

的像素点置为0，其它的像素点置为1。此时BWEDR_t中每个连通域只保存了为与连通域较中心的部分。将此结果再以边长为5的方形结构元素进行膨胀，得到最后的自适应腐蚀膨胀结果LastBW_t。利用连通域标记，得到它的连通域个数RegionNum_t。

(8)在水稻田中相机场景的四个角做标记，如图10所示，将四边形的观测场景划分为两个三角形，三角形1和三角形2，分别量取两个三角形的边长既可计算两个三角形的面积为S₁、S₂。设相机场景面积为SceneArea，于是SceneArea＝S₁+S₂一般每个稻穴中植入水稻数量为b，于是对于分割后的图像BW_t检测到的一亩植入水稻数量为：

${\mathrm{RiceNum}}_t={\mathrm{RegionNum}}_t\&times;b\&times;\frac{1}{0.0015\&times;\mathrm{SceneArea}}.$

本实例中，在检测阶段(4)、(6)、(7)的判别式中多次使用数组的均值与中值。这样可以进一步增强算法的自适应性，减少前面常数参数带来的影响。

本实施例检测阶段(4)中“最优膨胀结构元素要保证将腐蚀段的连通域合并在一起，而且保证不能过膨胀。随着结构元素的增大，检测到的连通域异常数量变化会开始的时候急剧下降，而后缓慢下降，急剧下降部分为合并大量腐蚀断的连通域造成的，缓慢下降部分为过膨胀时造成，于是两个阶段的转折点对应于最优膨胀结构元素的大小”，使用diff＞-15来检测转折点，这是本实例实验后的判别式。在应用中，不局限上述方法，只要检测到此曲线的转折点即可。

Claims (3)

1.一种水稻移栽后密度自动检测方法，具体为：

对水稻图像进行二值分割；

初始化R个半径不同的圆形状结构元素StrEroded；

利用R个圆形状结构元素StrEroded对分割得到的二值图像进行腐蚀得到R个腐蚀结果bwerode_tR；

创建L个边长不同的方形形状结构元素StrDilate；

利用L个方形形状结构元素StrDilate对腐蚀结果bwerode_tR进行膨胀，得到R×L个膨胀结果bwdilate_tRL；

检测各膨胀结果bwdilate_tRL的异常连通域数量；

选取异常连通域数量从急剧下降到缓慢下降的转折点对应的方形形状结构元素作为最优结构膨胀元素；

利用最优膨胀元素对腐蚀结果bwerode_tR进行膨胀；

从利用最优膨胀结构元素对腐蚀结果bwerode_tR进行膨胀的结果中选取连通域数量最多对应的圆形状结构元素为最优腐蚀结构元素；

利用最优腐蚀结构元素和最优膨胀结构元素对二值图像进行先腐蚀再膨胀得到BWED_t；

统计BWED_t的连通域个数RegionNum_t；

计算水稻密度 ${\mathrm{RiceNum}}_t={\mathrm{RegionNum}}_t\&times;b\&times;\frac{1}{0.0015\&times;\mathrm{SceneArea}},$ b为每个稻穴中植入水稻数量，SceneArea为相机场景面积。

2.根据权利要求1所述的水稻移栽后密度自动检测方法，其特征在于，所述

异常连通域是指满足公式(1)或(2)的连通域；

Distance≤K₁×MinDist    (1)

Distance为本连通域与其它连通域质心的最短欧式距离，MinDist为所有连通域的Distance的均值，0.3≤K₁≤0.6；

$\left\{\begin{array}{c}\mathrm{Pixnum}\&le;K_2\&times;\mathrm{MinPixnum}\\ \mathrm{Pixnum}\&GreaterEqual;K_3\&times;\mathrm{MinPixnum}\end{array}\right.---\left(2\right)$

Pixnum为本连通域的像素数，MinPixnum为所有连通域的Pixnum的均值，0.15≤K₂≤1，2≤K₃≤3。

3.根据权利要求1所述的水稻移栽后密度自动检测方法，其特征在于，在计算水稻密度步骤之前对所述BWED_t作去噪或/和连通域断开操作。

Images