一种机载LIDAR数据估算森林蓄积参数的方法

许 粲,江腾达

(1.安徽广播电视大学 六安分校,安徽 六安 237000；2.北科天绘科技有限公司,安徽 合肥 230000)

森林资源调查的主要目的是通过调查各类森林资源的分布和质量变化情况来进行有效的森林管理和决策分析。而干形数据是影响木材产量和质量的主要因子,传统的方法是通过人工对标准样地进行大量的实地勘测,并利用样本数据中单株树木的胸径、冠层直径和树高等参数来估算整个森林的林分信息,耗费大量的人力资源和时间成本,特别在一些人类难以深入的森林腹地,根本无法获取树木的准确数据。因此,随着数据更新的时效性越来越高,传统的作业模式已无法满足森林资源信息高效率、高精度的获取。在国外,无人机载激光雷达数据以其便捷、快速、低成本、并且可以获取高精度的树木高程等优势，被用于森林结构参数的估计中。陆地激光雷达扫描仪(TLS)可以获取高精度的树高、胸径等信息,但其获取的数据范围有限,不适合大面积勘测。机载激光雷达扫描仪(ALS)能够对大范围区域进行数据采集[1],在林业测量领域,它可以采集到树木的三维结构数据,对数据使用算法处理分析后,就能够建立树木的冠层结构模型,获得树高、树冠等特征参数[2],估算树胸径和冠层直径等信息。在单株树提取算法中,采用的数据源和方法多样,如利用波形数据[3]；方法上根据自底向上的分割策略,利用林分信息拓扑结构和权重估计2个树冠是否应归类为1株树来提取单独的树冠[4]；利用分水岭分割算法和分类结果相结合的方法进行单株树检测[5]等。

单株树木提取算法的准确性、点云密度大小以及森林分布和树冠结构等因素直接影响单株树木的提取精度。在采用图像分割方法之前,需要将三维点云数据转换成二维栅格数据,栅格图像数据中存在很多噪声点和间隙点,在此基础上通过窗口滑动检测局部最大值提取单株树木的冠点,并且将检测的树冠点作为标记图像用于分水岭分割算法时极易造成错分、多分问题。本文利用机载激光雷达采集的林区点云数据,通过数据预处理和图像分割等方法,提取试验区域内所有的树木和地面信息,在单株树的检测上,通过噪点消除法方法减少错分率。

1 算法原理

为验证机载激光雷达数据在林业调查应用中的可行性,设计方法流程如图1所示。

图1 方法流程图

1) 样地区域点云数据采集。数据采集需根据当地气候、高程起伏以及郁闭度条件,规划航线和航速,选择最佳飞行时间。利用无人机搭载激光雷达进行数据采集。

2) 预处理。将点云数据经纬度坐标进行UTM投影,转成以米为单位；利用Lastools(https://rapidlasso.com/lastools/)提供的Lasground工具,设置参数“-replace_z”来消除地形,获取归一化的数字高程模型nDSM[5]。也可以通过DSM(数字表面模型,表示树冠的顶部信息)与DTM(数字地形模型,表示地形数据)之间的差值来获取nDSM,nDSM又称为冠层高度模型CHM。

3) 种子点提取。利用局部最大值和噪点消除提取树木冠顶作为分水岭分割算法的标记图像，如图2(a)所示,点云转化成栅格数据时会出现较多的噪声点和间隙点,从而导致检测出的局部最大值不准确,将检测出的种子点作为标记图像(图2(c))进行分水岭分割时，会在局部小范围内出现多分现象(图2(e))。因此,本文将点云转化成栅格数据时的赋值进行调整,对所有落入固定栅格区域内的点进行比较,取最大值为当前灰度值(图2(b))。

RasterΩ(x,y)=max(Ω(x,y))

(1)

式中:RasterΩ(x,y)表示新生成的栅格图像,下标Ω表示局部格网区域；Ω(x,y)表示局部格网大小内落入的点云高程值。根据式(1)生成的标记图像以及在此基础上进行分水岭分割的结果图像分别如图2(d)、图2(f)所示,该方法有效地减少了多分、错分问题。

4) 分水岭分割算法提取单株树木。利用优化的标记图像进行分水岭分割运算,区分树的类别。再根据分割出的各个单株树木区域,反转换到三维点云空间。

5) 根据样地实测数据的林分情况选择胸径、树高、冠层直径之间的关系模型,计算森林蓄积量。

图2 标记点改进前后对比图

2 实例分析

2.1 样地区域

试验区域位于福建省三明市尤溪县,森林覆盖率高,分布密集。本次无人机载激光雷达共采集4组数据,分别如图3所示。样地大小为25.82m×25.82m(1亩),分别在每组数据的中间位置截取宽高相等的4块数据对各参数进行均值处理。试验区域占地总面积约为26.67hm2(400亩),采集数据耗费1天时间,相较人工测量(大于5 d)具有更高的作业效率。

2.2 数据处理

图3(c)数据3为例,分别给出原始点云数据(高程渲染)、nDSM、单株树提取结果,如图4所示。

图3 4组实验数据

图4 单株树检测中的过程数据

2.3 森林参数估计

利用激光雷达点云数据估算森林蓄积量,是在树检测算法完成后,已经识别出单株树的前提下,根据每株树的提取结果,定义每株树所在范围内的最大高程值为树高,冠径大小通过对位置范围内的所有点进行圆的拟合来估算。而胸径和蓄积量等参数则需要通过建立树高(H)、胸径(D)、冠径(C)三者之间的模型,根据特定树种、特定区域的特定公式来进行估算。根据相关学者的实验结果,认为大部分树种由树高和冠径单独估算胸径的精度没有结合树高和冠径来估算胸径的精度高。

由于本次实验未选取标准样地,没有进行树木相关参数的实地测量,因此无法拟合模型,这里采用地处经度25.11°,纬度61.19°区域内的针叶林拟合模型[6]:

(2)

式中:D,H,C分别表示胸径、树高和冠径。

这里,森林蓄积量采用平均实验形数法[7]进行估算,平均实验形数法的基本公式:

V=G1.3(H+3)fε

(3)

式中:V为森林蓄积量(m3)；G1.3为树干胸高断面积(m2)；H为树高(m)；fε表示平均实验形数。这里不考虑混合林等的影响,假设实验区的树种均为针叶林,因此取值为0.39。胸高断面积可表示为:

(4)

式中:a,b表示胸高处相互垂直的直径大小。当a=b时,用D表示,单位为m。

所以整片区域中总的森林蓄积量用下式表示(i表示第i株树):

(5)

2.4 结果与分析

分别对4组数据进行预处理,利用局部最大值和噪声消除提取标记图像,利用分水岭分割算法提取出单株树、对检测出的单株树进行圆拟合估算冠径、根据关系模型计算胸径以及森林蓄积量。其中,点云数据转成二维栅格图像时,设置的网格间距为0.2,局部最大值检测时,树最小半径设为1.5m；由于采集数据时周边区域的点很少,单个点也可能被误认为是树,因此,这里设置有效点数为20,即认为单株检测范围内的点数低于20个时,视为无效点。

表1给出每组数据的平均株数、平均树高、平均冠径、平均胸径以及平均森林蓄积量,并且与3组样地的部分人工实测数据进行比较。估算结果与部分样地的树高较为吻合,树高和冠径的大小直接影响平均胸径的结果,且一般情况下胸径是随着树高和冠径的增加而增加。树木的株数受局部最大值的检测窗口影响较大,因此,搜索窗口当前点与周边区域进行对比的横向宽度约为1～3m；特别在点云数据密度低、树种的冠形结构特征不明显,导致树木难以检测。

选取每组数据中的1个截取区域,表示单株树的检测结果和局部最大值点的对比情况。其中,绿色代表拟合圆的圆心和边界,红色表示局部最大值点,即为标记点,n表示检测的树木株数。如图5所示,点云数据质量好的地方除少量点稍有偏差,其余的局部最大值点均在拟合圆范围之内。而受点云密度和冠层结构影响,拟合的结果存在很大的偏差。

表1 各组数据森林参数信息

图5 树检测结果中拟合圆和局部最大值的分布图

3 结论

1) 传统人工测量工作受人力、地形和时间等因素制约,导致数据获取过程缓慢、精度低,而无人机载激光雷达可以实现快速、高精度地获取地面三维数据,可以为森林参数的估算提供数据支撑。

2) 单株树的检测结果受点云密度、冠层结构以及单株树提取算法的影响较大。特别当点云密度低、冠层形状散乱难以区分时,检测算法经常是失效的；树检测算法中,局部窗口选择过小会导致多分现象；过大,则出现少分现象。所以,一般情况下根据树的基本冠层宽度在[1,3]区间范围内进行设置。同时，本文通过噪点消除方法也进一步地减少了错分率。

3) 由于标准样地的实测数据不全,且缺乏完全对应的实验数据,无法建立树高-冠层直径-胸径之间的拟合模型,仅用1组已有的模型进行近似估算,而且没有与实测数据进行对比,无法给出算法的检测精度,后续应用中将对同步获取的实测数据进行对比分析。