基于无人机遥感的伐区造林坑穴参数提取方法

本发明涉及造林工程监测领域，具体涉及一种基于无人机遥感的伐区造林坑穴参数提取方法。

背景技术：

监测和评价造林工程的质量，是精准提升森林质量和实现森林效益的关键环节。林业生产方式落后、经营管理粗放、质量效益低下是我国林业存在的突出问题。只有转变林业经营方式，从传统粗犷林业向精准林业方向转换，才能有效提高森林质量，全面提升森林的多种功能，满足社会的多样化需求。

无人机遥感技术在植被信息提取和林业资源调查方面已经进行了大量的研究，各种方法被广泛应用于植被高度、植被覆盖度、造林规划设计、林分参数提取等方面。造林挖穴是植树造林的第一步也是造林过程中极其重要的一环，坑穴的数量和质量关乎植树造林的成活率，影响最终的造林成果质量。当前对造林坑穴的参数提取和质量评价主要依赖工作人员现场抽检，一方面工作量大效率低，另一方面客观准确性难以得到保证。随着图像处理技术的发展，基于无人机影像的人工智能和模板匹配技术在目标探测和图像自动解译等方面得到广泛应用。大量文献已经证明无人机遥感在各类目标探测和识别方面具有较好的实用价值。造林坑穴目标具有特定的形状、光谱和空间分布特征，相互之间具有较大的相似度。霍夫变换在提取特定形状目标方面有较多的应用，在伐区造林坑穴参数提取中具有应用潜力。

技术实现要素：

有鉴于此，本发明的目的在于提供一种基于无人机遥感的伐区造林坑穴参数提取方法，为造林坑穴质量评价提供自动化的解决方案，从而进一步拓展无人机遥感在造林质量检测中的应用潜力和价值。

为实现上述目的，本发明采用如下技术方案：

一种基于无人机遥感的伐区造林坑穴参数提取方法，包括以下步骤:

步骤s1:获取伐区造林坑穴的航拍影像数据，经处理得到dsm图像；

步骤s2:基于canny边缘检测方法，对dsm图像进行边缘检测，生成坑穴边界图像；

步骤s3:对坑穴边界图像圆形霍夫变换方法，得到圆形坑穴对象，将提取的坑穴对象个数作为坑穴数量提取结果；

步骤s4:根据坑穴数量提取结果，将圆形对象的直径作为坑穴的宽度；

步骤s5:根据步骤s3得到的坑穴提取结果和步骤s1生成的dsm，分别统计每个坑穴对象内部像素点高程值的最大值和最小值，两者的差值作为坑穴的深度。

进一步的，所述步骤s1具体为：采用无人机获取伐区造林坑穴的航拍影像数据，并对得到的航拍相片，使用摄影测量处理软件进行数据预处理，加入控制点并进行手工刺点，全自动空三加密生成影像密集匹配点云，生成实验区数字表面模型。

进一步的，所述无人机飞行设计航线为航向重叠率不低于90％，旁向重叠率不低于70％，空间分辨率为0.01-0.02m。

进一步的，所述步骤s2具体为：

步骤s21:对坑穴dsm图像进行高斯平滑滤波，计算公式如下：

式中,σ为高斯滤波器参数,它控制着平滑程度；x和y为图像中某一像素的位置；

步骤s22：对平滑后的图像i(x,y)计算梯度幅值和方向，计算公式如下：

θ(x,y)＝arctan(py(x,y)/px(x,y))

式中,x和y方向偏导数分别为:

px(x,y)＝(i(x,y+1)-i(x,y)+i(x+1,y+1)-i(x+1,y))/2

py(x,y)＝(i(x,y)-i(x+1,y)+i(x,y+1)-i(x+1,y+1))/2

步骤s23：遍历梯度幅值图m,在以m(x,y)为中心点的8邻域内判断,如果m(x,y)比梯度方向上相邻两个点的幅值都大，则令m1(x,y)＝m(x,y)，否则，令m1(x,y)＝0，m1为非极大值抑制结果；

步骤s24：分别使用高阈值thmax和低阈值thmin对图像m1进行判断,若m1(x,y)>thmax，判定该点是边缘点，若m1(x,y)<thmin,则判定该点不是边缘点，若thmax>m1(x,y)>thmin，则判定该点疑似边缘点，再进一步判断，若疑似边缘点的邻接像素中有边缘点，则认为该点也是边缘点，否则，认为该点为非边缘点。

进一步的，所述步骤s3中圆的表达式为：(x-a)²+(y-b)²＝r²，参数空间的表达式为：a＝x-rcosθ，b＝y-rsinθ，霍夫变换把求解图像空间中圆的曲线问题转换成求解参数空间中像素点最多的(a,b,r)参数对问题。其中，a和b表示圆心，r表示半径，具体包括以下步骤：

步骤s31设置最小半径rmin、最大半径rmax、半径步长rs、角度θ∈(0～2π)、角度步长θr，半径搜索步数cr＝(rmax-rmin)/rs，角度搜索步数为cθ＝2π/θr，初始化参数空间h(m,n,cr)，m为图像横坐标像素个数，n为纵坐标像素个数，h初始为全0矩阵；

步骤s32：边缘检测后的有效像素点坐标xi(xi,yi)，按照rs和θr，根据公式a＝x-rcosθ，b＝y-rsinθ进行网格搜索计算，得到ai，bi，ri，参数h(ai,bi,ri)累积加1，可以得到参数空间cr×cθ个h(ai,bi,ri)参数对；

步骤s33：重复步骤s32，计算图像空间中所有的边缘检测后的有效像素点，相同的h(ai,bi,ri)参数对进行加1；

步骤s34：计算完所有的有效像素点后，对h按照累加值大小排序。根据参数空间中h(ai,bi,ri)，确定圆心、半径，将圆输出到图像中，得到坑穴提取结果。

进一步的，所述坑穴深度为坑面点与坑底点高程之差，以圆形霍夫变换得到的圆心、直径分别作为窗口中心和窗口大小在原始dsm上进行局部像素点高程的极值探测，得到坑穴对象内部像素点高程值的最大值和最小值，两者的差值即为坑穴的深度参数。

进一步的，所述坑穴的深度参数提取公式如下：

式中,d为坑穴的深度，z(x,y)表示该坐标像素高程值，pv表示坑穴对象内包含的所有像素点。

本发明与现有技术相比具有以下有益效果：

本发明能够快速有效的获取伐区造林坑穴参数，为造林坑穴参数提取及质量检测提供一种新的方法。相比于现场人员抽查检测更加客观、科学和高效，是造林信息提取的最佳解决方案，为造林坑穴参数提取及质量检测提供一种新的方法，拓展了无人机遥感技术在造林质量评价方面的应用。

附图说明

图1为本发明实施例的方法流程示意图。

图2为本发明实施例的研究样地。

图3为本发明实施例的基于dsm局部圆形霍夫变换提取坑穴过程。其中a为dsm源数据，b为canny边缘检测结果。c为圆形霍夫变换探测结果，d为坑穴提取结果与dom(数字正射影像)叠加显示。

图4为本发明实施例的单个坑穴影像点云示例。

图5为本发明实施例的三种方法提取结果；其中a为圆形霍夫变换全区域结果，b为圆形霍夫变换局部结果，c为模板匹配局部结果，d为目视解译局部结果。

图6为本发明实施例的伐区造林坑穴实测和提取宽度散点图。

图7为本发明实施例的伐区造林坑穴实测和提取深度散点图。

具体实施方式

下面结合附图及实施例对本发明做进一步说明。

请参照图1，本发明提供一种基于无人机遥感的伐区造林坑穴参数提取方法，具体包括以下步骤：

步骤s1：利用无人机获取伐区造林坑穴的航拍影像数据，经处理得到数字表面模型(digitalsurfacemodel,dsm)；

步骤s2：基于canny边缘检测方法，对步骤s1生成的造林坑穴区dsm图像进行边缘检测，生成坑穴边界图像；

步骤s3：对步骤s2生成的坑穴边界图像应用圆形霍夫变换方法，得到圆形坑穴对象，将提取的坑穴对象个数作为坑穴数量提取结果；

步骤s4：基于步骤s3，将圆形对象的直径作为坑穴的宽度，即坑穴宽度等于圆形霍夫变换得到的圆形直径

步骤s5：基于步骤s3得到的坑穴提取结果和步骤s1生成的dsm，分别统计每个坑穴对象内部像素点高程值的最大值(max)和最小值(min)，两者的差值作为坑穴的深度。

在本实施例中，无人机飞行设计航线为航向重叠率不低于90％，旁向重叠率不低于70％，空间分辨率为0.01-0.02m；对得到的航拍相片，使用摄影测量处理软件pix4d进行数据预处理，加入控制点并进行手工刺点，全自动空三加密生成影像密集匹配点云，最后生成实验区数字表面模型(digitalsurfacemodel,dsm)。

在本实施例中，步骤s2具体包括以下步骤：

步骤s21：对坑穴dsm图像进行高斯平滑滤波，计算公式如下：

式中,σ为高斯滤波器参数,它控制着平滑程度；x和y为图像中某一像素的位置。

步骤s22：对平滑后的图像i(x,y)计算梯度幅值和方向，计算公式如下：

θ(x,y)＝arctan(py(x,y)/px(x,y))

式中,x和y方向偏导数分别为:

px(x,y)＝(i(x,y+1)-i(x,y)+i(x+1,y+1)-i(x+1,y))/2

py(x,y)＝(i(x,y)-i(x+1,y)+i(x,y+1)-i(x+1,y+1))/2

步骤s23：遍历梯度幅值图m,在以m(x,y)为中心点的8邻域内判断,如果m(x,y)比梯度方向上相邻两个点的幅值都大，则令m1(x,y)＝m(x,y)，否则，令m1(x,y)＝0，m1为非极大值抑制结果；

步骤s24：分别使用高阈值thmax和低阈值thmin对图像m1进行判断,若m1(x,y)>thmax，判定该点是边缘点，若m1(x,y)<thmin,则判定该点不是边缘点，若thmax>m1(x,y)>thmin，则判定该点疑似边缘点，再进一步判断，若疑似边缘点的邻接像素中有边缘点，则认为该点也是边缘点，否则，认为该点为非边缘点。

在本实施例中，步骤s3中圆的表达式为：(x-a)²+(y-b)²＝r²，参数空间的表达式为：a＝x-rcosθ，b＝y-rsinθ，霍夫变换把求解图像空间中圆的曲线问题转换成求解参数空间中像素点最多的(a,b,r)参数对问题。其中，a和b表示圆心，r表示半径。具体包括以下步骤：

步骤s31：初始化。设置最小半径rmin、最大半径rmax、半径步长rs、角度θ∈(0～2π)、角度步长θr，半径搜索步数cr＝(rmax-rmin)/rs，角度搜索步数为cθ＝2π/θr，初始化参数空间h(m,n,cr)，m为图像横坐标像素个数，n为纵坐标像素个数，h初始为全0矩阵；

优选的，在本实施例中rmin＝17，rmax＝30，rs＝1，θr＝0.1；

步骤s32：霍夫变化求解参数对。边缘检测后的有效像素点坐标xi(xi,yi)，按照rs和θr，根据公式a＝x-rcosθ，b＝y-rsinθ进行网格搜索计算，得到ai，bi，ri，参数h(ai,bi,ri)累积加1，可以得到参数空间cr×cθ个h(ai,bi,ri)参数对；

步骤s33：计算参数对累积值。重复步骤s32，计算图像空间中所有的边缘检测后的有效像素点，相同的h(ai,bi,ri)参数对进行加1；

步骤s34：输出识别图像。计算完所有的有效像素点后，对h按照累加值大小排序。根据参数空间中h(ai,bi,ri)，确定圆心、半径，将圆输出到图像中，得到坑穴提取结果。

在本实施例中，步骤s5中坑穴深度为坑面点与坑底点高程之差。以圆形霍夫变换得到的圆心、直径分别作为窗口中心和窗口大小在原始dsm上进行局部像素点高程的极值探测，得到坑穴对象内部像素点高程值的最大值和最小值，两者的差值即为坑穴的深度参数。提取公式如下：

式中,d为坑穴的深度，z(x,y)表示该坐标像素高程值，pv表示坑穴对象内包含的所有像素点。

实施例1：

本实施例选择的造林坑穴样地位于福建省将乐县某伐区(26°51′13″n，117°34′50″n)，该林场位于福建省西北部，地处武夷山脉东南部、金溪河畔，多为中、低山地貌，以低山丘陵为主，海拔高平均在400-800m，最高处海拔为1403m，而最低海拔为140m。研究区整体坡度约为25°，面积约3712m²，图2为研究区示意图。

实验数据使用大疆精灵4pro型无人机进行遥感影像数据的采集，该型无人机搭载1英寸2000万像素cmos，等效35mm焦距，携带方便，在复杂环境下也可快速执行影像获取任务。因测区地形起伏大且无人机不具备仿地形飞行的能力，结合实际将测区划分为两部分设计航线来获取数据，其目的是防止出现地面分辨率变化过大导致图像拼接效果差的情况。设计航线为航向重叠率90％，旁向重叠率70％，飞行相对高度约为36m，空间分辨率为0.01m。为了保证成果的平面和高程精度以及将成果统一至cgcs2000(chinageodeticcoordinatesystem2000)坐标系下,现场采用华测i70型gnss接收机高精度获取11个地面控制点(groundcontrolpoint,gcp)，rtk(real-timekinematic)模式下的平面精度为±(8+1×10^-6×d)mm，高程精度为±(15+1×10^-6×d)mm。

使用主流无人机摄影测量处理软件pix4d进行数据预处理，得到实验区正射影像(digitalorthophotomap,dom)、数字表面模型(digitalsurfacemodel,dsm)及无人机摄影测量点云产品。对使用gcp和不使用gcp生产的dsm分别进行精度验证，得到的标准误差(rootmeansquareerror,rmse)分别为0.056m与0.158m。根据林业生产要求，标准坑穴深度0.400m，误差小于15％，使用gcp生产的dsm满足量测的精度要求。

图3中的a是数据预处理生成的局部dsm图像，经过canny边缘检测处理，坑穴边界像素大于0，其他像素值全部归零，从而得到准确的坑穴边界信息，如图3中b所示；基于canny算子得到的边界图像应用圆形霍夫变换找圆，将物理空间长度参数米转换为图像空间中像素数，得到坑穴提取结果，如图3中c所示；将坑穴提取结果与dom叠加显示，如图3中d所示。提取结果的数量即为坑穴数量参数，提取圆形直径即为坑穴的宽度参数。

由图4单个坑穴的三维影像点云可知，坑穴深度为坑面点与坑底点高程之差。分别统计每个坑穴对象内部像素点高程值的最大值和最小值，两者的差值即为坑穴的深度参数。

选取坑穴局部提取结果如图5，图5中d为目视解译结果，将其作为正确的识别结果，对比分析其他两种提取效果。可以看出，基于dom模板匹配法(如图5中c)相比于基于dsm圆形霍夫变换法(如图5中b)有较多的坑穴误提和漏提现象；基于dsm圆形霍夫变换法则与目视解译结果最为接近。

利用检验指标，分别计算两种方法提取识别精度，如表1所示。结果表明，基于dsm圆形霍夫变换坑穴提取法精度最高，其正确率为95.15％，误检率为4.85％，漏检率为6.56％。该法相较于基于dom模板匹配法提取的结果，正确率提升明显，漏检率得到进一步的控制，更加符合质检需要。

实际人工采集20个坑穴的宽、深参数对基于dsm圆形霍夫变换坑穴宽、深提取结果进行精度验证，由实测宽度与提取宽度、实测深度与提取深度线性拟合结果如图6、图7分析，二者的r²分别达到0.93和0.92，rmse分别为1.02cm和1.67cm，提取参数和实测值之间呈现明显的相关性，提取结果较为可靠。

对比三种方法在提取坑穴数量和宽深参数的耗时、精度、工作量以评估处理速度和适用性，由表2中可以看出，相比于人工目视解译方法的耗时，其余两种方法在处理速度上显然更具优势；在考量精度、工作量和自动化程度基础上，基于dsm圆形霍夫变换法既可提取数量又可提取宽度和深度信息，因而具有最高的适用性。

注：√表示可实现；×表示不可以实现

对于山区地形来说，因为地形起伏、太阳光线变化以及周边高大树木遮挡的影响，最终生产得到的dom数据可能会存在色彩不均匀的现象。基于模板匹配法提取坑穴对象要求dom色彩分布均匀，当坑穴目标位于阴影下或者曝光过度区域则容易导致漏检和误检现象的产生，从而使得整体漏检率和误检率偏高。此外，模板质量对模板匹配法提取坑穴的精度影响较大，样本的选取要求较高，这在一定程度上降低了处理效率和精度。

本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等)上实施的计算机程序产品的形式。

本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上所述，仅是本发明的较佳实施例而已，并非是对本发明作其它形式的限制，任何熟悉本专业的技术人员可能利用上述揭示的技术内容加以变更或改型为等同变化的等效实施例。但是凡是未脱离本发明技术方案内容,依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与改型，仍属于本发明技术方案的保护范围。