结合三维密集点云的无人机影像大豆覆盖度提取

何海清，严椰丽,2，凌梦云，杨勤锐，陈 婷，李 麟

（1. 东华理工大学测绘工程学院，南昌330013；2. 宁波市阿拉图数字科技有限公司，宁波 315000；3. 中国地震局第二监测中心，西安 710054；4. 东华理工大学水资源与环境工程学院，南昌330013；5. 江西裕丰智能农业科技有限公司，赣州341309）

0 引 言

作物冠层覆盖度是反映植物生长状况的常用参数，它直接决定作物对光的截流量，进而影响光合作用效率，是对作物长势进行分析的常用指标之一。快速、准确地获取大豆冠层覆盖度对于大豆农田精细化管理和产量估测等具有重要价值。国内外研究人员对植被覆盖度提取进行了广泛而深入的研究，早期研究仅限于利用手持摄像机采集单幅影像来提取作物覆盖度，不适用于大范围作物覆盖度提取。近年来，以无人机为平台的低空遥感飞速发展，使大面积田间影像数据的获取更加便捷和多样化。相比传统的影像获取手段，无人机数字摄影测量具有成本低、操作灵活等优点，可获取高分辨率的遥感影像并能解算出高精度的三维点云数据，能为大范围农情信息监测提供了新的技术手段。因此，利用无人机遥感技术进行植被覆盖度提取已成为研究热点之一，且取得了良好的应用效果，具有代表性的研究方法主要分为3类；一是仅基于可见光植被指数进行二维遥感影像分割，实现植被与背景的分割以提取农作物覆盖度；二是结合可见光植被指数与其他分类方法如：面向对象分类法、Otsu自动阈值分割法、监督分类法、像元二分法等；三是通过颜色空间变换后再分割作物与土壤，实现农作物信息提取。然而，前期研究表明，利用常用的归一化的绿红差值指数、可见光大气阻抗植被指数、绿叶指数、比率指数等植被指数难以通过设定通用的阈值来最优化提取植被信息。此外，以上研究主要采用二维影像来提取植被覆盖度，受影像分辨率、光照条件、背景复杂度等诸多因素制约，导致在包含低矮杂草、地膜等复杂背景下难以精确地提取作物覆盖度。

针对当前仅利用二维遥感影像难以从光照不均匀和复杂背景中准确获得大豆冠层覆盖度的问题，本研究结合无人机摄影测量三维点云进行大豆覆盖度提取，利用改进的可见光植被指数来准确提取植被信息，并采用经过局部阈值算法分割后的数字表面模型来进行空间分析，以剔除低矮杂草等噪声影响，进而实现三维空间内大豆覆盖度提取。花芽分化期是大豆生长较关键的时期，此生长期大豆尚未完全封垄、覆盖度相对较低、田间杂草生长旺盛，获取该时期的大豆覆盖度对于作物长势分析、产量估测等具有十分重要的现实意义。因此，本研究选取该时期大豆作为研究对象进行作物冠层覆盖度提取。

1 材料与方法

1.1 研究区概况

为验证本研究方法在起伏地形下大豆种植覆盖度提取的有效性，选取中部地区包含平地和地形起伏的具有代表性大豆种植区进行试验。研究区域为湖北省京山市永隆镇（30°44′41″N，112°50′5″E）的大豆种植区，区位图见图1a，该地区总面积约100 km，地势北高南低，属暖温带东亚半湿润季风区，季风气候，四季分明，光照充足，雨量充沛，年平均降水量1 179 mm，年内各个季节气候差异明显，年平均气温16.3 ℃，主要种植大豆和小麦等粮食作物。

为评估本研究方法在不同时期、不同杂草混杂程度、地形起伏下剔除低矮杂草等复杂背景的有效性，在研究区域内选取具有代表性的3个地块进行试验，该区域数字正射影像图（Digital Orthophoto Map，DOM）和三维密集点云高程渲染图见图1b～1d：地块1位于研究区南部地形较为平坦区域，面积约为1 076 m，田间存在较多细小、低矮杂草等复杂背景的大豆花芽分化初期；地块2位于研究区南部地形较为平坦区域，面积约1 325 m，大豆处于花芽分化中期，包含更多、更密的低矮杂草；地块3位于研究区北部，面积约3 378 m，存在明显的地形起伏变化，田间存在较多、较密的混杂低矮杂草等复杂背景，大豆生长期为花芽分化中期。采用等距穴播栽培法进行种植，植株行距和穴距分别约为35、20 cm，同一时期大豆长势情况大致相同。

1.2 遥感影像获取与预处理

本研究试验数据采集时间为2020年6月4日、14日（晴朗、云少、风速小于10 m/s），通过小型四旋翼无人机平台（DJI御Mavic Air，大疆，中国）拍摄获取的立体影像，该无人机飞行平台主要参数如表1所示。本次试验无人机飞行高度设定为10 m，地面分辨率约为0.4 cm，飞行速度为1.2 m/s，镜头垂直向下拍摄，航向和旁向重叠度均设置为80%，每幅影像的分辨率像素为4 056×3 040，影像以24位真彩色.jpge格式存储。

表1 无人机飞行平台主要参数 Table 1 Parameters of UAV flying platform

由于低成本无人机携带的相机几何结构不稳定，导致其拍摄的影像存在较大的畸变。首先，基于相机参数对影像进行畸变差校正，并重采样生成无畸变的影像；其次，利用运动恢复结构（Structure from Motion，SfM）算法与区域网平差恢复影像摄站姿态、重构地物稀疏三维点云；最后，利用半全局匹配（Semi-Global Matching，SGM）算法生成研究区包含真彩色纹理信息的稠密三维点云，并通过微分纠正法制作正射影像。本研究试验得到的DOM和三维密集点云渲染图见图1b～1d所示。

图1 研究区位置及正射影像和密集点云 Fig.1 Location, DOM and dense point cloud of three fields in the study area

1.3 研究方法

本研究的技术路线如图2所示，主要包括以下5个步骤：1）通过摄影测量技术如相机标定、SfM、空三、SGM密集匹配等完成2个试验区无人机影像空中三角测量，并生成高分辨率正射影像和附带真彩色纹理信息的三维密集点云；2）采用伽马增强的绿叶指数（Gamma-enhanced Green Leaf Index，GGLI）精确提取植被信息；3）采用最佳结构元的局部阈值分割法分割摄影测量数字表面模型；4）基于植被信息与数字表面模型分割结果的相交分析来剔除低矮杂草等复杂背景噪声；5） 通过目视判读构建混淆矩阵来评价本研究大豆覆盖度提取精度。

图2 技术路线图 Fig.2 Workflow of the proposed method

在本研究中，为了增强植被指数特征并抑制不相关的背景噪声，通过增强植被光谱反射率来构建伽马增强绿叶指数GGLI以达到精确提取植被信息的目的。GGLI的定义如式（1）所示。

式中GGLI映射为8 bit影像灰度值，取值范围为[0, 255]；表示伽马增强系数，经多次试验以达到最佳植被提取效果确定系数值为2.5；、、分别为红、绿、蓝真彩色的3个分量值。在本研究中，将GGLI0> 的部分视为植被，其他部分视为土壤背景，能较完整地提取绿色植被信息。然而，仅依赖二维影像提取植被信息虽能有效剔除非绿色植被背景，但难以剔除低矮杂草等复杂背景噪声。

为有效剔除低矮杂草影响，本研究提出一种结合三维信息的局部阈值影像分割方法来精确提取大豆覆盖度。首先，利用三维密集点云生成研究区不规则三角网（Triangulated Irregular Network，TIN）和数字地表栅格图，则作物与背景物体的高度区别可映射为栅格图灰度差异。由于农田地形起伏和作物高度差异，导致生成的地表栅格图灰度分布不均匀，难以通过设定全局阈值来精准分割目标物。因此，本研究采用一种基于局部阈值的影像分割方法来实现大豆作物分离，该方法通过顶帽变换和Otsu算法获得局部阈值函数对影像进行二值化，能有效克服灰度分布不均匀的影响，准确分离作物与背景。

顶帽变换定义为从原影像中减去其数学形态学开运算后的影像，能有效解决影像中背景纹理与光照不均匀、目标物不突出等问题。顶帽变换的数学表达式如式（2）所示。

式中img、(img)分别为原影像和顶帽变换后的影像；img°表示结构元对影像img的开运算操作。

式中img⊗ 和img⊕ 分别表示当结构元的原点位于像素点坐标(,)时，用结构元在该位置进行腐蚀和膨胀操作；和为结构元的相关参数。

通过开运算操作可剔除比结构元更小的噪声点，而保留面积较大的区域。在开运算中，假设有一个球形结构元，推动该结构元沿影像三维灰度曲面的下侧面来回移动，直至移动位置覆盖整个下侧面，此时球体任何部分达到的最高点即可构成开运算之后的曲面。

然后，利用Otsu算法对顶帽变换结果进行二值化处理，该算法具有直观性和易实现等优点。该算法利用影像灰度特性，可将影像分割成背景和感兴趣区2个部分，该方法又被称为最大类间方差法，即通过Otsu算法分割后的背景与感兴趣区域的类间方差达到最大。

假设灰度级的 Otsu算法分割阈值为() = （0 ＜＜-1），将影像分割结果分为两类：背景（）和目标（）。、分别表示影像灰度值在范围[ 0,]和[+ 1,- 1]内的所有像素，则每一类出现的概率（）、（），以及各类的平均灰度值()、()如式（7）～式（10）所示。 同时，由灰度级0至累加均值 ()及整幅影像的平均灰度计算如式（11）和式（12）所示。

基于上述计算结果，阈值最大化的类间方差σ()和全局方差σ如式（13）和式（14）所示。

最佳阈值的最大类间方差σ()数学表达式如式（15）所示。

在确定最佳阈值后，将输入影像二值化为背景和目标区域（即0和1），其数学表达式如式（16）所示。

式中=0,1,2,… ,-1；=0,1,2,…,-1。

由于仅利用二维可见光谱影像提取大豆冠层覆盖度无法剔除低矮杂草等非大豆植被，本研究结合摄影测量三维密集点云和相交分析过滤低矮杂草等噪声。将提取后的植被影像与局部阈值分割影像作为对象进行相交分析，则相交部分被视为精化后的大豆区提取结果，并计算大豆覆盖度。该相交分析数学表达式如式（17）所示。

式中为公共像素点集合，(,)、 (,)分别为提取后的植被影像与局部阈值分割后的影像在像素点位置(,)的数值。∩表示影像的相交运算。若和在同一像素位置都为非零值，则将该像素点位置逻辑值设置为“是”，反之，将该像素点设置为“否”。

作物冠层覆盖度（Canopy Coverage，CC）通常定义为统计范围内作物地上部分（包括叶、茎、枝）的垂直投影面积所占的百分比，其计算如式（18）所示。

式中为影像中大豆冠层像素点个数总和；为影像像素点个数总和，即×。

1.4 精度评价

为评价本研究方法的大豆覆盖度提取精度，采用相对误差（Relative Error，RE）衡量大豆覆盖度提取结果与大豆覆盖度参考值之间偏差程度，反映提取方法的可靠性，RE越小，大豆覆盖度提取精度越高，计算公式如下：

式中、C分别表示大豆覆盖度提取结果与大豆覆盖度参考值。

同时，采用目视判读构建混淆矩阵的方法来计算总体分类精度（Overall Accuracy，OA）和Kappa系数评价。

总体分类精度（Overall Accuracy，OA）定义为被正确分类的特征点数量占总特征点数量的百分比，其数学表达式如式（20）所示。

式中、分别为作物和背景被正确分类的样本数；为总样本数。

Kappa系数是一个用于检验分类结果与实际结果一致性的指标，可用于衡量分类的效果，其取值范围为[-1，1]，通常大于0，其值越接近1，说明分类效果越好。其数学表达式如式（21）和式（22）所示。

式中、分别为作物和背景被误分的样本数。

2 结果与分析

本研究结合可见光谱与三维密集点云进行大豆覆盖度提取试验，根据豆苗形态特征，经多次试验，最终采用半径为120个像素点的圆形结构元对研究区数字表面模型栅格影像进行开运算操作，可较好地从复杂背景中分离出大豆，该方法不易受地形起伏影响，并能有效去除低矮杂草等噪声。最后，结合可见光谱和三维信息对提取后的植被影像进行相交分析，实现大豆覆盖度精确提取。为便于对比分析，利用ArcGIS矢量化工具，通过人工判读手动提取大豆覆盖度作为参考值，地块1、2、3的大豆覆盖度参考值分别为42.51%、36.88%、33.63%。研究区3个地块大豆覆盖度提取结果见图3，白色、黑色区域分别表示大豆和背景信息。

由图3可知，基于二维影像GGLI的方法可较好地提取无杂草的大豆信息。通过公式（18）计算，地块1、2、3的大豆覆盖度分别为46.33%、41.75%、38.81%，但无法准确提取包含杂草背景的大豆覆盖区域，致使提取结果存在较多的噪声，3个地块提取大豆覆盖度相对误差RE分别为8.99%、13.20%、15.40%，可见随着地块1、2、3的杂草混杂程度与地形起伏等背景复杂性加剧，基于二维影像GGLI的方法提取大豆覆盖度误差也逐步变大。而本文研究方法能把大部分低矮杂草剔除，地块1、2、3的大豆覆盖度分别为41.16%、37.65%、34.17%，其大豆覆盖度相对误差RE分别为3.17%、2.09%、1.61%，与人工判读的参考值较为吻合，可见本文研究方法提取大豆覆盖度并未受到杂草混杂程度与地形起伏的影响。由于本研究方法受限于影像像素光谱数值，相比人工勾画的大豆覆盖区域，以像素为单元的提取方法存在提取大豆区域不连片、面积细小等现象。

图3 试验地块大豆覆盖度提取结果对比 Fig.3 Comparisons of soybean coverage extraction results

此外，采用目视判读构建混淆矩阵并计算分类精度来定量化评价本研究方法针对大豆覆盖度提取的效果。基于前期研究经验，针对每一试验地块，随机选取1 000个样本点，采用等分原则平均分配样本点数量，即500个作物点，500个背景点。首先，在正射影像中分别手动标记作物点和背景点，再叠加相交分析结果进行目视判读，并构建混淆矩阵，见表2。由目视判读可知，在选取的500个作物点中，在地块1、2、3中，被错误分类为背景点数分别为8、6、9；在选取的500个背景点中，在地块1、2、3中，被错误分类为作物点数分别为6、5、7，可见本文研究方法得到错误的作物点和背景点数较少。由混淆矩阵计算得到大豆覆盖度提取的总体精度（OA）和Kappa系数见表3，由定性分析和前述相对误差指标与混淆矩阵定量化精度评价结果可知，相比仅依赖二维影像的提取方法，本研究方法可有效地从无人机低空摄影测量数据中分离大豆与背景信息，能得到较高的分类精度和实现大豆覆盖度精准提取。

表2 试验地块大豆覆盖度提取结果混淆矩阵 Table 2 Confusion matrix of soybean coverage extraction results in the study fields

表3 试验地块大豆覆盖度提取结果总体精度和Kappa系数 Table 3 Overall accuracy and Kappa coefficient of soybean coverage extraction results in the study fields

为进一步验证本研究方法的有效性和适用性，选用几种具有代表性的基于二维图像分割的覆盖度提取方法来进行对比分析，包括支持向量机（Support Vector Machine，SVM）分类、结合Lab颜色空间变换与Kmeans分割法、双峰阈值法分割、基于过绿指数（Excess Green Index，ExG）、植被提取颜色指数（Color Index of Vegetation Extraction，CIVE）和红绿比指数（Green Red Ratio Index，RGRI）分割，分割结果见图4，经SVM、Lab-Kmeans、双峰阈值、ExG、CIVE和RGRI二维图像分割后，田垄上低矮杂草噪声显然未能与大豆作物分离，无法得到准确的大豆覆盖度。相比而言，本文结合三维信息的局部阈值分割方法能较好地剔除低矮杂草等复杂背景噪声，以精准提取大豆覆盖度。

图4 试验地块大豆覆盖度不同提取方法结果对比 Fig.4 Extraction results of different extraction methods for soybean coverage in study fields

这几种方法所提取的大豆覆盖度对比见图5、总体精度OA（精度评价方法同上）见图6，仅基于二维图像分割方法（M1～M6）提取大豆覆盖度超过参考覆盖度5%，而本研究提出基于伽马增强的绿叶指数分割方法M7在3个地块的大豆覆盖度提取结果与覆盖度参考值相比误差均小于3%，总体精度OA要优于方法M2～M6和略低于M1。由图5对比3个地块大豆覆盖度参考值可知，相对于地块1，地块2大豆生长更为旺盛、地块3地形起伏变化，方法M1～M6提取的大豆覆盖度受不同生长期与地形起伏影响较大和参考值差距更大。在图6中，方法M1-M2的总体精度OA相对M3～M6略高1%，其原因在于大豆覆盖度较大致使被误分的作物样本点相对少一点，但仍然与目视判读参考值差异较大。而本研究同时结合可见光谱和三维信息提取大豆覆盖度与覆盖度参考值最为接近，误差小于1.5%，相比而言，3个试验地块大豆覆盖度提取总体精度OA值最高（优于98%）。基于视觉定性和指标定量分析，试验验证了本研究方法结合二维影像和三维点云能较好地剔除细小、低矮杂草等复杂背景噪声和准确地提取出花芽分化期大豆冠层覆盖度，对于不同时期、不同杂草混杂程度、不同地形起伏背景下大豆覆盖度提取具有较强的鲁棒性。此外，在同一硬件条件下，各方法耗时见表4，可见SVM耗时最长，相比而言，基于指数的大豆覆盖度提取效率远高于基于分割的方法（高68%以上），而本研究方法大豆覆盖度提取耗时与基于指数的方法较为接近。

表4 不同方法耗时对比 Table 4 Time consuming comparison of different methods s

图5 地块1、2、3不同方法大豆覆盖度提取统计 Fig.5 Statistics of bean coverage extraction by different methods in field 1, 2, 3

图6 试验地块不同方法大豆覆盖度提取的总体精度比较 Fig.6 Comparisons of overall accuracy of bean coverage extraction by different methods in study fields

3 结 论

本研究针对仅依赖二维遥感影像提取大豆覆盖度难以剔除杂草等复杂背景干扰的问题，提出了一种结合三维点云进行大豆覆盖度提取的方法。以花芽分化期大豆的可见光无人机影像为研究对象，采用伽马增强植被指数提取植被信息，结合最佳结构元的局部阈值分割法二值化附带三维信息的数字表面模型，并通过相交分析来过滤非大豆植被等低矮地物，实现不同时期、不同杂草混杂程度、不同地形起伏背景下大豆覆盖度精确提取。主要结论如下：

1）基于伽马增强的可见光植被指数（Gamma- enhanced Green Leaf Index，GGLI）可增强植被指数特征和抑制不相关的背景噪声，达到有效剔除非植被信息的目的，从而大幅提高绿色植被提取精度。

2）本研究采用最佳结构元的局部阈值二值化附带三维信息的数字表面模型能有效剔除低矮杂草等复杂背景噪声，对于光照和地形起伏变化具有较强的鲁棒性。

3）结合植被提取结果与三维点云信息相交分析能进一步剔除非大豆复杂背景干扰，实现大豆覆盖度的精确提取，相关试验总体精度达到98%以上，相比支持向量机、结合Lab颜色空间变换与Kmeans分割法、双峰阈值法等常用方法效率高68%以上，在精度和效率方面明显优于仅利用二维影像的覆盖度提取方法。