CN112669333A - 一种单木信息提取方法 - Google Patents

Abstract

本发明公开了一种单木信息提取方法，包括以下步骤：S1、扫描待测森林区域，得到原始点云；S2、基于统计分析技术，对原始点云进行去噪处理，得到去噪点云；S3、对去噪点云内的采样点进行滤波式分类，分类出地面点和非地面点；S4、根据非地面点构建待测森林区域数字表面模型DSM，采用TIN法根据地面点构建待测森林区域数字高程模型DEM，并计算数字表面模型DSM和数字高程模型DEM两者间的差值，得到冠层高度模型CHM；S5、基于分水岭算法，分割得到单木的树高和冠幅信息；本发明解决了传统技术对于单木信息的提取存在计算方法不足、难以克服地形劣势、结果不精确的缺点，从而导致过分割和欠分割的问题。

Description

一种单木信息提取方法

技术领域

本发明涉及激光雷达点云数据处理领域，具体涉及一种单木信息提取方法。

背景技术

随着遥感传感器技术及航空和航天平台技术的迅速发展，尤其是LiDAR技术的出现，为森林资源调查带来了新活力，使得林业资源调查工作更加的现代化，综合化，精细化。LiDAR主动向探测目标发射高频率的激光脉冲，并直接获取地物表面的精准三维空间坐标和回波信息等，在估测森林冠层高度及林木冠层空间结构方面具有独特优势，提高森林关键结构参数特别是冠层垂直结构参数的估测精度。然而，随着森林资源调查精细化程度的不断提高以及林业应用的不断深入，现有的一些技术已经不能满足相关工作的需求。

发明内容

针对现有技术中的上述不足，本发明提供的一种单木信息提取方法解决了传统技术对于单木信息的提取存在计算方法不足、难以克服地形劣势、结果不精确的缺点，从而导致过分割、欠分割的问题，所以提出了一种基于单木信息提取的方法来更好的克服地形等劣势条件，让单木信息的提取结果变的更准确。

为了达到上述发明目的，本发明采用的技术方案为：一种单木信息提取方法，包括以下步骤：

S1、通过LiDAR激光雷达，扫描待测森林区域，得到原始点云；

S2、基于统计分析技术，对原始点云进行去噪处理，得到去噪点云；

S3、通过渐进加密三角网滤波算法，对去噪点云内的采样点进行滤波式分类，分类出地面点和非地面点；

S4、根据非地面点构建待测森林区域数字表面模型DSM，采用TIN法根据地面点构建待测森林区域数字高程模型DEM，并计算数字表面模型DSM和数字高程模型DEM两者间的差值，得到冠层高度模型CHM；

S5、基于分水岭算法，根据冠层高度模型CHM，分割得到单木的树高和冠幅信息。

进一步地，步骤S2包括以下分步骤：

S21、计算得到原始点云中每一个采样点与其最近邻域的n个采样点的距离值，并根据距离值，计算得到距离均值μ和距离标准差σ，其中n为正整数；

S22、遍历每一个采样点与其最近邻域的n个采样点的距离值，在遍历的每个环节中，删除当前采样点的最近邻域的n个采样点中距离值大于μ+std_mul×σ的采样点，得到去噪点云，其中std_mul为去噪阈值。

进一步地，步骤S3包括以下分步骤：

S301、以1米长度作为空间格网尺寸，对去噪点云的三维空间进行格网化划分；

S302、在格网化划分的三维空间中，选取每个格网内最低处的去噪点云采样点作为起始种子点，建立三角网；

S303、在去噪点云中随机选取一采样点记为被分析点；

S304、查询被分析点水平面投影所落入的三角网，记为被分析三角网；

S305、计算被分析点到被分析三角网的垂直距离d，以及被分析点到被分析三角网的三个顶点与被分析三角网所在平面所成角度的最大值α；

S306、判断d＜1.4且α＜8°是否成立，若是，则跳转至步骤S307，若否，则跳转至步骤S308；

S307、标记被分析点为地面点，并将其加入被分析三角网，跳转至步骤S308；

S308、判断去噪点云中是否存在未被选取作为被分析点的采样点，若是，则跳转至步骤S309，若否，则跳转至步骤S310；

S309、在去噪点云中未被选取作为被分析点的采样点中随机选取一采样点记为被分析点，并跳转至步骤S304；

S310、记去噪点云中所有未被标记为地面点的采样点为非地面点。

上述进一步方案的有益效果为：渐进加密三角网滤波算法对于不连续地形，尤其是植被密集的城市地区，能够得到最佳的滤波结果。

进一步地，步骤S4中构建待测森林区域数字高程模型DEM包括以下分步骤：

A1、对去噪点云的三维空间的水平面进行栅格划分，栅格边长为1m；

A2、以每个栅格中高程最小的地面点的值作为该栅格的栅格值；

A3、判断是否存在因内部无地面点而无法被赋值的栅格，若是，则跳转至步骤A4，若否，则得到数字高程模型DEM，结束；

A4、采用TIN法，根据已赋值的栅格，对因内部无地面点而无法被赋值的栅格进行插值，使得每个栅格均被赋值，得到数字高程模型DEM。

上述进一步方案的有益效果为：采用TIN法根据地面点构建待测森林区域数字高程模型DEM能够很好的表示复杂地形，避免了因内插而造成的精度损失。

进一步地，步骤S4中构建待测森林区域数字表面模型DSM包括以下分步骤：

B1、对去噪点云的三维空间的水平面进行栅格划分，栅格边长为1m；

B2、设置参考点(x₀，y₀)，对栅格(i，j)插入内插点(x₀+i，y₀+j)，其中，(i，j)为栅格的坐标，(x₀+i，y₀+j)为内插点的坐标；

B3、构建矩形检索函数，并基于矩形检索函数检索每一个内插点在矩形检索函数搜索区域内的非地面点，判断矩形检索函数搜索区域内非地面点是否达到额定数量，若是，则得到被分析内插值点p，若否，则调整矩形检索函数搜索区域的范围，直到矩形检索函数搜索区域内非地面点达到额定数量，得到被分析内插值点p；

B4、计算每个非地面点到被分析内插值点p的反距离权P_i，其中，

u为梯度参数，d_pi为被分析内插值点p与第i个非地面点之间的距离；

B5、根据反距离权P_i，计算被分析内插值点p的高程值Z_p，得到构建的待测森林区域数字表面模型DSM。

上述进一步方案的有益效果为：采用反距离加权法能够快速、高效的将非地面点利用起来生成DSM，以便于后续的计算。

进一步地，步骤B5中高程值Z_p的计算公式为：

其中，Z_i为被分析内插值点p的邻域内的第i个非地面点，n为邻域内的非地面点数量。

进一步地，步骤S5包括以下分步骤：

S51、对冠层高度模型CHM计算梯度值；

S52、根据梯度值，绘制梯度图；

S53、采用分水岭算法对梯度图进行第一次分水岭分割，得到第一次分割图像；

S54、对第一次分割图像的区域最小值进行定位，得到最小值邻近的区域；

S55、对最小值邻近的区域进行第二次分水岭分割，得到第二次分割图像；

S56、对第二次分割图像进行去除极端值处理，并提取单木信息，得到单木的树高和冠幅信息。

上述进一步方案的有益效果为：基于CHM影像梯度图加两次分水岭分割的算法能够在地形起伏和植被覆盖剧烈变化的区域，有效精确的分割出单木信息，且能够将错误分割，即过分割、欠分割现象的出现概率降到最低。

进一步地，步骤S51中计算梯度值的公式为：

其中，G为梯度值，G_x为经横向边缘检测的图像灰度值，G_y为经纵向边缘检测的图像灰度值。

综上，本发明的有益效果为：在林业领域的研究中，若能实现从海量机载LiDAR点云信息中分离出高精度的单株木的三维信息，可以有效地减少在传统森林调查中对时间、人力及经济成本的消耗与浪费。若能实现高精度的单木点云分割，则可以从分割结果中直接得到树高、冠幅等树木结构信息。此外，作为最重要的森林资源管理手段，在森林资源清查中应用LiDAR技术及相关单木分割方法，将大大节约森林资源清查的时间及人工成本。

附图说明

图1为一种单木信息提取方法的流程图；

图2为被分析三角网的示意图；

图3为测试区域的DEM、DSM和CHM在直方图均衡化后显示的图；

图4为梯度图；

图5为第一次分割效果图；

图6为最小值区域图；

图7为第二次分割效果图；

图8为第一张分割标记效果图；

图9为第二张分割标记效果图；

图10为单木因子图；

图11为单木位置、树高和冠幅图。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

如图1所示，一种单木信息提取方法，包括以下步骤：

S1、通过LiDAR激光雷达，扫描待测森林区域，得到原始点云；

S2、基于统计分析技术，对原始点云进行去噪处理，得到去噪点云；

林区的机载LiDAR点云数据不仅包含地物信息，而且还存在少量噪声点，需要对其进行去噪处理。常见的噪声包括高位粗差和低位粗差。通过选择合适的参数，可以移除噪点，提高数据质量。本次测试使用统计分析技术，从一个点云数据中集中移除测量噪声点。对每个测量噪声点的邻域进行统计分析，剔除不符合一定标准的邻域点。具体步骤如下：

步骤S2包括以下分步骤：

S21、计算得到原始点云中每一个采样点与其最近邻域的n个采样点的距离值，并根据距离值，计算得到距离均值μ和距离标准差σ，其中n为正整数，在本实施例中n为10；

S22、遍历每一个采样点与其最近邻域的n个采样点的距离值，在遍历的每个环节中，删除当前采样点的最近邻域的n个采样点中距离值大于μ+std_mul×σ的采样点，得到去噪点云，其中std_mul为去噪阈值，在本实施例中，std_mul为5。

S3、通过渐进加密三角网滤波算法，对去噪点云内的采样点进行滤波式分类，分类出地面点和非地面点；

LiDAR点云数据滤波在这里是指在进行了去噪后的点云数据中分离出地面点和非地面点的过程，目的是为了获取DEM和DSM。本实施例区域位于茂密林区，点云数据中主要为森林冠层激光点、地面点，所以滤波在这里就是为了分离森林冠层激光点和地面点。针对去噪后的LiDAR点云数据，本实施例采用的是改进的渐进加密三角网滤波算法分类地面点。

步骤S3包括以下分步骤：

S301、以1米长度作为空间格网尺寸，对去噪点云的三维空间进行格网化划分；

S302、在格网化划分的三维空间中，选取每个格网内最低处的去噪点云采样点作为起始种子点，建立三角网；

S303、在去噪点云中随机选取一采样点记为被分析点；

S304、查询被分析点水平面投影所落入的三角网，记为被分析三角网；

S305、计算被分析点到被分析三角网的垂直距离d，如图2所示，以及被分析点到被分析三角网的三个顶点与被分析三角网所在平面所成角度的最大值α；

S306、判断d＜1.4且α＜8°是否成立，若是，则跳转至步骤S307，若否，则跳转至步骤S308；

S307、标记被分析点为地面点，并将其加入被分析三角网，跳转至步骤S308；

S308、判断去噪点云中是否存在未被选取作为被分析点的采样点，若是，则跳转至步骤S309，若否，则跳转至步骤S310；

S309、在去噪点云中未被选取作为被分析点的采样点中随机选取一采样点记为被分析点，并跳转至步骤S304；

S310、记去噪点云中所有未被标记为地面点的采样点为非地面点。

S4、根据非地面点构建待测森林区域数字表面模型DSM，采用TIN法根据地面点构建待测森林区域数字高程模型DEM，并计算数字表面模型DSM和数字高程模型DEM两者间的差值，得到冠层高度模型CHM，如图3所示；

相比于传统的遥感光学影像数据，LiDAR系统发射的激光脉冲能够穿透林区树冠直接获取林区地面点。利用LiDAR估测森林参数和进行森林资源调查时，首要的基础工作就是获取林区地面激光点云和冠层激光点云，进而分别建立林区数字高程模型(DEM)和数字表面模型(DSM)。DSM表征林区冠层形态，即由滤波后的非地面点构建的模型。DEM是数字地面模型的一个分支，能够较为真实地模拟林下地形，即根据滤波过后地面点进行构建的地表模型。CHM是冠层高度模型，是DEM和DSM的差值。

步骤S4中构建待测森林区域数字高程模型DEM包括以下分步骤：

A1、对去噪点云的三维空间的水平面进行栅格划分，栅格边长为1m；

A2、以每个栅格中高程最小的地面点的值作为该栅格的栅格值；

A3、判断是否存在因内部无地面点而无法被赋值的栅格，若是，则跳转至步骤A4，若否，则得到数字高程模型DEM，结束；

A4、采用TIN法，根据已赋值的栅格，对因内部无地面点而无法被赋值的栅格进行插值，使得每个栅格均被赋值，得到数字高程模型DEM。

步骤S4中构建待测森林区域数字表面模型DSM包括以下分步骤：

B1、对去噪点云的三维空间的水平面进行栅格划分，栅格边长为1m；

B2、设置参考点(x₀，y₀)，对栅格(i，j)插入内插点(x₀+i，y₀+j)，其中，(i，j)为栅格的坐标，(x₀+i，y₀+j)为内插点的坐标；

B3、构建矩形检索函数，并基于矩形检索函数检索每一个内插点在矩形检索函数搜索区域内的非地面点，判断矩形检索函数搜索区域内非地面点是否达到额定数量，若是，则得到被分析内插值点p，若否，则调整矩形检索函数搜索区域的范围，直到矩形检索函数搜索区域内非地面点达到额定数量，得到被分析内插值点p；

B4、计算每个非地面点到被分析内插值点p的反距离权P_i，其中，

u为梯度参数，d_pi为被分析内插值点p与第i个非地面点之间的距离；

式中指数u＞0，而且实验证明，当u＞2时，导致曲面在数据点附近相对比较平直，而在两个数据点之间一个很小的区域内有很大的梯度，当u＜2时，导致曲面相对平缓，没有起伏，当u＝2时，不但容易计算，也比较符合实际地形变化规律，因此，实际应用中常常取u＝2。

B5、根据反距离权P_i，计算被分析内插值点p的高程值Z_p，得到构建的待测森林区域数字表面模型DSM。

步骤B5中高程值Z_p的计算公式为：

其中，Z_i为被分析内插值点p的邻域内的第i个非地面点，n为邻域内的非地面点数量。

S5、基于分水岭算法，根据冠层高度模型CHM，分割得到单木的树高和冠幅信息。

步骤S5包括以下分步骤：

S51、对冠层高度模型CHM计算梯度值；

步骤S51中计算梯度值的公式为：

其中，G为梯度值，G_x为经横向边缘检测的图像灰度值，G_y为经纵向边缘检测的图像灰度值。

S52、根据梯度值，绘制梯度图，如图4所示；

S53、采用分水岭算法对梯度图进行第一次分水岭分割，得到第一次分割图像，如图5所示；

S54、对第一次分割图像的区域最小值进行定位，得到最小值邻近的区域，如图6所示，左图为最小值位置图，右图为邻近图；

S55、对最小值邻近的区域进行第二次分水岭分割，得到第二次分割图像，如图7所示；

图8为分割标记效果图，图中白圈黑点为盆地的中心，黑线为分水岭。

S56、对第二次分割图像进行去除极端值处理，并提取单木信息，得到单木的树高和冠幅信息。

除极端值处理，进行连通区域面积统计，如果大于设置的最大连通区域250或小于最小连通区域6，则舍弃。最后的分割结果如下图9所示。

提取单木信息。在软件arcmap中，依据分割线，叠加CHM，提取单木的树高；按照分割区域，求取冠幅，这里求的是几何重心到边缘线的平面距离。最终结果下图10、11所示。

Claims (8)

1.一种单木信息提取方法，其特征在于，包括以下步骤：

S1、通过LiDAR激光雷达，扫描待测森林区域，得到原始点云；

S2、基于统计分析技术，对原始点云进行去噪处理，得到去噪点云；

S3、通过渐进加密三角网滤波算法，对去噪点云内的采样点进行滤波式分类，分类出地面点和非地面点；

S4、根据非地面点构建待测森林区域数字表面模型DSM，采用TIN法根据地面点构建待测森林区域数字高程模型DEM，并计算数字表面模型DSM和数字高程模型DEM两者间的差值，得到冠层高度模型CHM；

S5、基于分水岭算法，根据冠层高度模型CHM，分割得到单木的树高和冠幅信息。

2.根据权利要求1所述的单木信息提取方法，其特征在于，所述步骤S2包括以下分步骤：

S21、计算得到原始点云中每一个采样点与其最近邻域的n个采样点的距离值，并根据距离值，计算得到距离均值μ和距离标准差σ，其中n为正整数；

S22、遍历每一个采样点与其最近邻域的n个采样点的距离值，在遍历的每个环节中，删除当前采样点的最近邻域的n个采样点中距离值大于μ+std_mul×σ的采样点，得到去噪点云，其中std_mul为去噪阈值。

3.根据权利要求1所述的单木信息提取方法，其特征在于，所述步骤S3包括以下分步骤：

S301、以1米长度作为空间格网尺寸，对去噪点云的三维空间进行格网化划分；

S302、在格网化划分的三维空间中，选取每个格网内最低处的去噪点云采样点作为起始种子点，建立三角网；

S303、在去噪点云中随机选取一采样点记为被分析点；

S304、查询被分析点水平面投影所落入的三角网，记为被分析三角网；

S305、计算被分析点到被分析三角网的垂直距离d，以及被分析点到被分析三角网的三个顶点与被分析三角网所在平面所成角度的最大值α；

S306、判断d＜1.4且α＜8°是否成立，若是，则跳转至步骤S307，若否，则跳转至步骤S308；

S307、标记被分析点为地面点，并将其加入被分析三角网，跳转至步骤S308；

S308、判断去噪点云中是否存在未被选取作为被分析点的采样点，若是，则跳转至步骤S309，若否，则跳转至步骤S310；

S309、在去噪点云中未被选取作为被分析点的采样点中随机选取一采样点记为被分析点，并跳转至步骤S304；

S310、记去噪点云中所有未被标记为地面点的采样点为非地面点。

4.根据权利要求1所述的单木信息提取方法，其特征在于，所述步骤S4中构建待测森林区域数字高程模型DEM包括以下分步骤：

A1、对去噪点云的三维空间的水平面进行栅格划分，栅格边长为1m；

A2、以每个栅格中高程最小的地面点的值作为该栅格的栅格值；

A3、判断是否存在因内部无地面点而无法被赋值的栅格，若是，则跳转至步骤A4，若否，则得到数字高程模型DEM，结束分步骤；

A4、采用TIN法，根据已赋值的栅格，对因内部无地面点而无法被赋值的栅格进行插值，使得每个栅格均被赋值，得到数字高程模型DEM。

5.根据权利要求1所述的单木信息提取方法，其特征在于，所述步骤S4中构建待测森林区域数字表面模型DSM包括以下分步骤：

B1、对去噪点云的三维空间的水平面进行栅格划分，栅格边长为1m；

B2、设置参考点(x₀，y₀)，对栅格(i，j)插入内插点(x₀+i，y₀+j)，其中，(i，j)为栅格的坐标，(x₀+i，y₀+j)为内插点的坐标；

B3、构建矩形检索函数，并基于矩形检索函数检索每一个内插点在矩形检索函数搜索区域内的非地面点，判断矩形检索函数搜索区域内非地面点是否达到额定数量，若是，则得到被分析内插值点p，若否，则调整矩形检索函数搜索区域的范围，直到矩形检索函数搜索区域内非地面点达到额定数量，得到被分析内插值点p；

B4、计算每个非地面点到被分析内插值点p的反距离权P_i，其中，

u为梯度参数，d_pi为被分析内插值点p与第i个非地面点之间的距离；

B5、根据反距离权P_i，计算被分析内插值点p的高程值Z_p，得到构建的待测森林区域数字表面模型DSM。

6.根据权利要求5所述的单木信息提取方法，其特征在于，所述步骤B5中高程值Z_p的计算公式为：

其中，Z_i为被分析内插值点p的邻域内的第i个非地面点，n为邻域内的非地面点数量。

7.根据权利要求1所述的单木信息提取方法，其特征在于，所述步骤S5包括以下分步骤：

S51、对冠层高度模型CHM计算梯度值；

S52、根据梯度值，绘制梯度图；

S53、采用分水岭算法对梯度图进行第一次分水岭分割，得到第一次分割图像；

S54、对第一次分割图像的区域最小值进行定位，得到最小值邻近的区域；

S55、对最小值邻近的区域进行第二次分水岭分割，得到第二次分割图像；

S56、对第二次分割图像进行去除极端值处理，并提取单木信息，得到单木的树高和冠幅信息。

8.根据权利要求7所述的单木信息提取方法，其特征在于，所述步骤S51中计算梯度值的公式为：

其中，G为梯度值，G_x为经横向边缘检测的图像灰度值，G_y为经纵向边缘检测的图像灰度值。

Images