基于地基激光雷达的活立木材积提取算法

熊妮娜，王佳

(1.北京市城市管理研究院，北京 100028；2.北京林业大学精准林业北京市重点实验室，北京 100083；3.北京林业大学测绘与3S技术中心，北京 100083)

传统的立木材积主要依靠胸径尺、测高器等工具对标准样地进行调查，以标准木胸径、树高查询材积表进行估测，该方法效率低、人为干扰因素多，受材积表精度影响较大[1]。目前全站仪作为立木材积精准量测的仪器已经开展了相关实验和应用研究，该方法技术成熟，但同样存在操作复杂、不能海量获取坐标点等问题。而地基激光雷达技术作为现代测量的新兴技术手段之一，近年来被广泛地应用在单木计测和森林调查各个方面，激光雷达扫描可以数秒内一次性采集到上万点，在获取目标三维坐标位置方面具有先天优势。目前在林业中应用于森林生态系统、森林资源调查的各领域，包括树种、树高、胸径、生物量等单木属性信息的提取，叶面积指数的估算，冠层空隙度，冠层辐射，冠层结构，叶面积分布，树干曲线，竞争指数等[2]。随着地基激光雷达体积逐渐减小，更加便于携带，同时生产成本的降低带来价格的降低，地基激光雷达将会更加普遍地应用于林业调查中。

目前国内外学者开展了利用地基激光雷达提取单木参数的一系列研究。Thies等[3]、Astrup等[4]提出了一种基于立木的地面激光扫描数据重建树干三维表面的方法和算法。Tansey等[5]根据地面激光扫描仪数据估算了科西嘉岛松林地的树木材积。邓向瑞等[6]利用地基激光雷达对活立木进行了精确扫描测量，利用软件提取树木的直径、树高、冠幅等测树因子，利用材积模型计算立木材积，并与伐倒木的材积作对比。王剑等[7]通过激光点云数据构建一系列圆柱模型片段叠加而成树干的三维模型，该方法对树干的体积计算也提供了一种便利。蒋佳文等[8]采用地面激光雷达进行多站点扫描获取立木的点云信息，提取有关点云分布的特征参数，建立基于特征参数的材积模型。参考国内外学者的研究，本研究基于地基激光雷达点云数据，改进Hough变换方法，通过提取不同高度处的树干直径，运用断面积区分求积法计算树干材积，以期为无损精准测木提供新的方法。

1 数据处理与获取

1.1 外业数据采集

本研究采用FARO 120型地基激光雷达，扫描仪稳定可靠，测量距离为0.6～120 m，分辨率为0.1 mm，数据获取率≤50 800像素/s，25 m内误差≤2 mm，水平和垂直视野范围为360°和320°。

在树干胸高处朝南或朝北方向贴标靶纸，选取样木周围无遮挡的地方均匀放置3个参考球，选择能同时见到3个参考球的均匀分布的3个位置设站，理想的测站位置间隔角度为120°(以待扫描树为参考)，安装地基激光雷达并连接好相关设备，对扫描范围、扫描分辨率等参数进行设置后，进行树木扫描。本研究设置扫描区域为水平方向360°，垂直方向155°，扫描分辨率为10 000点/周，每株样木需从不同角度扫描至少3次，耗时约10 min/株。

1.2 内业数据处理

本研究点云数据处理参考张浩[9]、张冬等[10]采用的枝叶分离方法，主要包括4个步骤，结果如图1所示：

1)点云数据预处理。将地基激光雷达扫描数据导入计算机中，通过其配套软件Faro Scene，对点云数据进行加载、配准、滤波、剔除、导出等内业处理。打开Faro Scene软件，直接把测站数据加载到软件中，第一次打开软件时，扫描数据的Mio扫描点数值默认为62，这个数值是根据计算机的内存自动显示的，本研究涉及的点云提取使用的计算机内存为2 G，经过在不同内存的计算机中进行点云数据提取实验，结果显示，默认的Mio扫描点的数值均能满足精度要求。

2)树木点云数据主方向计算。树木点云数据主方向计算需要借助近邻点坐标信息及近邻点相应的法向量信息，原始的树木点云只包含点云坐标信息，本研究借助于主成分分析进行点云法向量计算；法截线拟合法计算主方向对噪声点具有一定的免疫力，是目前计算点云主方向较好的一种方法，计算出的枝干部分的点云主方向与枝干中轴线平行，冠层部分的点云主方向杂乱无章。

图1 内业处理过程Fig. 1 The interior work process

3)树木点云数据枝叶分割。首先根据点云数据各点及相应的主方向构建圆柱，圆柱半径定为2倍采样间距，圆柱高度定为10倍采样间距。然后将圆柱包围盒内的点投影到圆柱轴线上，根据圆柱内的点与圆柱轴线的夹角、投影点在轴线上的分布密度进行判断，区分冠层叶子点和枝干点。研究中分割时阈值∂1设为0.7，阈值∂2设为30°，阈值∂3设为0.8，其结果较好。

4)枝干骨架特征点提取。枝干部分依据点云测地图划分水平集，基于最小二乘法对水平集进行二维圆拟合，同时通过目视判断对拟合结果进行修正。

2 立木材积提取算法

2.1 任意处直径提取算法

本研究利用自动检测任意处直径的方法，首先截取任意处直径的点云数据水平切片，并进行栅格化，然后再使用圆拟合结合Hough变换，这是一个特征提取技术在特定形状不完善的情况下由一个决策程序找到对象的规则形状的方法，最终将圆检测出来，进而计算直径。该方法分为两个步骤：

图2 胸高部分水平切片截取算法和胸径量测算法流程Fig. 2 The flow chart of DBH of individual trees retrieval algorithm

1)任意处直径部分水平切片截取方法及步骤

任意处直径的点云数据水平切片是从构成单木的点云数据中选取高于任意处直径0.1 m至低于直径处0.1 m的点云数据。以Matlab实现算法，具体过程如图2所示。单木树干任意直径处上下0.1 m的小切片被截取下来，将此部分点云数据投影到xy平面，以便进行下一步边缘检测，并确定检测圆的半径。

2)任意直径处量测算法及步骤

在林业测量中，常遇到直径切面为不规则圆的情况，故本研究在不规则的直径切面点云投影数据上检测最能满足截面轮廓的标准圆，该检测圆的直径即视为树干任意处直径。圆的检测和半径的计算通过Hough变换方法，由于Hough变换检测圆程序需要读入图片，所以首先对水平切片投影后的数据进行栅格化处理，结合格网大小，换算出真实圆的直径[11-12]。变换圆的方程为x2+y2+2ax+2by+c=0，用两个累加器A(a，b)和B((a2+b2-c)1/2)来进行统计可能的圆心和半径。

本研究采用了一种改进的Hough变换算法，通过降低随机采样的点数，以及利用采样点的特征信息来判断是否进行累积，来避免大量的无目标的随机采样与无用累积[13-14]，使算法性能大大提高。具体算法如图2所示。经测试，将最小圆周上边缘点的个数定义为2 000较为合适。由于原始的点云数据是5个站点配准叠加的数据，数据量较大，如果参数设置太小，会导致检测圆的效果不佳，检测出多个小圆。第一次随机选取3个边缘点，将这3个边缘点代入变换圆公式x2+y2+2ax+2by+c=0，用两个累加器A(a，b)和B[(a2+b2-c)1/2]来进行统计可能的圆心和半径。然后清除属于该检测圆的所有边缘点。由于这里视树干直径为一个标准圆，故所有小圆的圆心位置视为树干直径标准大圆的边缘点，进行边缘检测，当某次循环结束但产生边缘点的统计个数达到或超过阈值的构造圆时，则圆检测结束。检测结果不是单个圆，而是以多个小圆圆心为边缘点组成多个半径一致、圆心位置相近的圆。所得结果对于获得树干直径无影响，若想获得圆心位置，将得到的多个圆圆心位置均值可看作所求圆心半径，即为单木的位置。最后结合栅格图像格网大小，得到单木的真实任意树干处直径。

以实验提取的一棵山杨为例，提取的树干不同高度的点云切片见图3a，选取胸径处的点云切片数据为例进行后面的直径提取见图3b，选取位置水平切片投影后的单木直径点云数据分布见图3c，栅格化后的水平切片见图3d，利用Hough变换检测出的圆效果见图3e，圆心坐标、半径长度和真实单木直径在Matlab软件中可自动计算出来。对本研究的单木，其圆半径为56，再结合预先设置的栅格格网大小0.1 cm，树干直径即为56×0.1×2=11.2 cm。

图3 利用Hough变换检测出圆的过程Fig. 3 Detect circle by using Hough transform

2.2 单木材积计算模型

本研究采用中央断面积区分求积法计算单木材积，将树干按一定长度(本研究采用2 m)分段，量出每段中央直径和最后不足一个区分段梢头底端直径，利用中央断面近似求积式，求算各分段的材积(V1，V2，V3，…，Vn)，并合计(V)。

V=V1+V2+V3+…+Vn+V′=g1l+g2l+g3l+

(1)

2.3 实验精度验证方法

对于研究立木材积精度验证最可靠的方法为：将点云数据提取的立木材积和树木伐倒后计算的材积相对比，但由于受限于政策原因，无法对树木进行此类破坏性试验。因此，本研究利用全站仪活立木材积测量方法对实验树木提取立木材积作为材积近似真值，全站仪无损测量立木的原理如下，首先用胸径尺测量立木树干的地径D0和胸径D1处(距离地面1.3 m处)的胸高h1，然后用全站仪测量仪器中心到胸径D1处的斜距S，再用全站仪将胸径以上的树干分成n段并测量相应位置的水平角αi和天顶距γi(下角标i为分段号)，后运用三角高程原理计算树高H，并按照圆台累加法计算得到立木树干的总材积V，具体方法如图4所示。目前已有研究表明，全站仪立木材积提取方法获取的材积与伐倒木真实材积误差在5%以内[15-16]，可以近似作为真值检验其他方法。

注：H为树高，m；h1, …, hi, …, hn为各分段高度，m；V为总材积，m3；V1, …, Vi, …, Vn为各分段材积，m3；D0、D1分别为地径和胸径，cm；D2, …, Di, …, Dn为各分段直径，cm；α1,…, αi, …, αn为各分段处的水平角，(°)；γ1, …, γi, …, γn为各分段处的天顶距，(°)；S为仪器中心到胸径处的斜距，m。 图4 全站仪测量立木材积原理Fig. 4 Schematic diagram for measuring standing tree volume using total station

e=V激光V全站

(2)

(3)

(4)

式中：V全站为全站仪提取材积数据，m3；V激光为地基激光雷达提取材积数据，m3。

3 结果与分析

本研究实测共56棵山杨树，采用现场直接量测的树木胸径做为参考值，将地基激光雷达扫描提取的胸径与其进行对比，其结果见表1。其中绝对误差最大的为0.023 3 m，最小值为0.008 0 m，没有绝对误差的有9棵，占总数的16.07%；绝对误差为正值的共有35棵，占总数的62.50%；为负值的为12棵，占总数的21.43%。从相对误差来看，最大的为8.82%，最小的为0.00%，平均相对误差为1.81%，可以看出地基激光雷达扫描提取胸径与现场直接量测的树木胸径十分接近，精度较高。

针对实测的56棵山杨树，运用本研究地基激光雷达扫描提取的材积与全站仪提取材积进行对比，其结果见表1。其中，绝对误差最大的为 0.070 6 m3，最小值为-0.110 2 m3，绝对误差为正值的共有37棵，占总数的66.1%，为负值的为19棵，占总数的33.9%。从相对误差来看，最大的为15.07%，最小的为0.06%，平均相对误差为5.86%，可以看出地基激光雷达扫描提取材积与全站仪提取材积十分接近，误差较小。同时分析实验结果发现，地基激光雷达扫描提取的材积常常会比全站仪提取的稍大，特别是个别树木出现了10%以上的相对误差。通过对误差较大树木的点云数据观察和分析发现，产生该问题的原因在于利用全站仪观测时，由于人为操作，如某一树高处，遮挡严重观测不到树干特征点，则可变换位置进行观测，能够比较好地提取干形，而地面激光雷达则只是固定的三站，扫描的点云数据无法完整包括不同高度处树干的特征信息，特别是树叶较为密集处提取直径时，可能将部分树叶和枝干的点云数据参与计算，进而造成了点云数据提取的材积普遍大于全站仪提取的材积。

表1 地基激光雷达扫描提取的胸径和材积Table 1 Data analysis of DBH and volume based on ground-based laser scanner

将样木实测数据导入SPSS软件中，建立两种方法的一元线性模型，公式如下：

V全站=1.02V激光-0.012

(5)

计算该线性模型的相关系数R为0.998，决定系数为0.996，说明两种方法具有显著的相关性。该模型回归系数的标准误差为0.022，说明模型回归系数较好，进一步验证了两种方法相关性极高。

由此，可以得出地基激光雷达扫描方法可以作为活立木材积提取方法。

4 结论与讨论

1)本研究采用改进的Hough变换方法提取树干任意处直径，避免了以前随机Hough 变换的无效累积，提高了算法性能。实验采用直接测量的胸径值作为参考，通过地基激光雷达量测的胸径值平均相对误差为1.81%；实验通过采用基于地基激光雷达提取材积与全站仪提取材积进行比较，材积十分接近，平均相对误差为5.86%，并且存在显著相关性。有效验证了基于地基激光雷达提取活立木材积的可行性。

2)本研究以地基激光雷达这一快速、准确获取被测物体三维空间结构信息的技术为基础，实现了对样木的三维重构。研究选取了56棵山杨树进行地基激光雷达扫描，样木之间遮挡不严重，枝叶密度不高，树干点云容易提取，扫描效果较好。未来的研究中应选择多个树种开展实验，验证该方法的可行性。

3)与全站仪获取立木材积方法比较，地基激光雷达虽然同样需要多次设站扫描树木，后期还需对点云数据各种处理、研建模型等过程，精度相似，过程更麻烦。但应用地基激光雷达技术的真正意义在于一次扫描可以提取包括活立木材积在内的诸多森林调查因子，例如叶面积指数、冠层空隙度、枝叶密度等全站仪难以获取的参数，而树木健康状况、树种识别等需要纹理信息的参数全站仪更是无法做到，在这些方面地基激光雷达的应用已经得到了很好的实验验证。同时，地基激光雷达比全站仪操作简单，测量时间短，全站仪通常测量1棵立木需要30 min以上，地基激光雷达单次扫描3～5 min，1棵树10 min即可完成，效率有很大的提升。而内业数据处理部分，虽然地基激光雷达相对麻烦，但是通过编程开发计算程序以后可以重复使用。因此，利用本研究提出的方法提取活立木材积，在森林资源监测中地基激光雷达有广阔的应用前景。