范永祥 冯仲科 申朝永 闫 飞 苏珏颖 王 蔚

(1.季华实验室， 佛山 528000； 2.北京林业大学精准林业北京市重点实验室， 北京 100083;3.贵州省第三测绘院， 贵阳 550004)

0 引言

森林资源及其监测所获取的可靠信息可以为不同层级、不同目标的林业部门制定可持续发展政策、投资决策等提供支持，从而为森林培育、森林采伐、气候变化影响评估、火灾模型和碳储量估算等奠定基础[1]。森林资源清查是获取森林资源信息的重要途径。由于不同层级的使用者对森林资源信息的使用目的不同，所关心的森林资源属性也不尽相同[2]。为满足清查目的，需要根据森林资源及所要求的精度制定不同的调查策略。样地调查是一种重要的森林资源清查方法，其通过在研究区域合理布设样地，然后对样地的森林属性进行调查、统计，从而反演出整个区域总体或平均的森林属性[3]。通常样地直接调查的属性包含树种、胸径、树高及立木位置等，生物量、碳储量等则可通过以上属性构建的反演模型进行估计[4]。传统样地调查中，使用一些简单的仪器进行森林调查(如使用胸径尺或卡尺完成胸径测量，使用罗盘及皮尺等完成立木位置测量)，这些方法不仅耗时耗力，且无法克服观测者主观因素的干扰。随着测距、测角等相关技术的发展，Laser-relascope[5]、电子测树枪[6]等设备提供了同时测量树木胸径、树高、位置等的一体化解决方案，但仍然无法克服安装复杂、观测难、受主观因素干扰等问题。

随着传感器技术的发展以及计算机运算能力的增强，通过点云提取样地单木及林分信息的方法得到快速发展。点云通常通过数字摄影测量技术、ToF(Time of flight)技术等获取。数字摄影测量系统通常以相机作为观测传感器，利用相机记录的纹理信息及角度重建被观测目标的三维采样及三维模型，但由于森林中植被遮挡严重，很难单纯利用摄影测量技术高效完成样地调查；此外，由于林下纹理过于复杂，很难形成样地完整三维点云[7-9]。ToF相机[10-13]由于其分辨率低、测距精度低且易受自然光线等影响，也不足以高效完成样地调查[14-15]。

相比于普通相机及ToF相机，激光雷达(Light detection and ranging，LiDAR)具有高效、不受自然光线影响等特性，是近年来比较流行的样地清查解决方案[16]。地基激光雷达(Terrestrial laser scanning，TLS)单次扫描由于仅有一个扫描视角，立木遮挡将无法获取到较大样地的完整三维数据，这导致所提取样地属性精度降低或无法获取到完整的样地属性[17-18]。多次扫描虽然避免了数据遗漏问题，但合并各站点扫描数据具有一定难度，且在大样地中布设太多站点耗时费力，TLS不适合在较大的样地中使用[19-20]。移动激光雷达系统(Mobile laser scanning，MLS)通常指以GNSS(Global navigation satellite system)、IMU(Inertial measurement unit)、激光雷达为主要传感器的运动平台，其通过在样地中动态采集样地数据，后处理时利用GNSS+IMU组合获取的姿轨曲线拼接激光雷达数据[21-22]；但通常森林郁闭度较高，受树冠遮挡影响，GNSS接收机很难在林下正常工作。即时定位及建图(Simultaneous localization and mapping，SLAM)技术的出现使得在林下无GNSS信号的区域定位成为可能[23-26]。随着SLAM算法的改进，单台多线激光雷达采集数据即可完成定位及建图，单激光雷达相比于TLS、MLS等具有设备简单、成本低等优势，其中一种典型的LiDAR SLAM算法为LOAM(LiDAR odometry and mapping in real-time)[27-30]，该算法以单帧数据中面特征、线特征为单帧去畸变、配准等过程的特征，具有计算速度快等优势；由于森林中线、面特征较少，单帧去畸变过程可能引入较大误差，导致配准精度低、构建点云噪声大等，故很难将该算法直接用于森林调查中。

针对以上问题，本文以LOAM为基础构建森林样地调查系统。在SLAM系统工作流中引入二次去畸变、二次配准等模块以提高去畸变、配准的鲁棒性及精度；使用32线激光雷达扫描4块32 m×32 m的森林样地，利用改进的LOAM系统完成样地建图；从点云提取立木位置及胸径，并通过参考数据及LOAM估计结果对比，完成LiDAR SLAM系统在森林样地中建图精度的间接评估。

1 LOAM算法及其变体

在处理连续单帧LiDAR数据时(以第k帧Pk为例)具体步骤为：

(1)基于曲率提取单帧点云数据中的线、面特征。

(2)假设Pk扫描期间雷达做匀速运动，即若tk时刻至tk+1时刻雷达位姿变换量为

(1)

则tk时刻至ti时刻雷达位姿变换量可线性内插为

(2)

图1 LOAM及其变体工作流程图Fig.1 Workflow of LOAM and its alternatives

从LOAM的结构中可以看出其主要结构为雷达里程计，并不包含传统SLAM系统的后端(回环检测及图优化)，故并不适用于较大的建图场景。

LeGO-LOAM(Lightweight and ground-optimized LiDAR odometry and mapping on variable terrain)作为一种重要的LOAM变更体[28]，主要进行了如下改进：①将单帧点云中地面及周围点云进行了分类。②在特征提取时，提取地面点云中的面特征、周围点云中的线特征。③在去畸变时，利用面特征完成竖直方向的线元素及翻滚角、俯仰角的优化，利用线特征完成其他元素的优化。④以ICP(Iterative closest point)算法为基础构建了回环检测功能，并添加了图优化功能。相比于LOAM，该系统在室内或有建筑物的室外等场景建图精度更佳，且适用于更大规模的建图场景。

2 LiDAR.SLAM算法改进

相比于室内或有建筑物的室外场景，森林中具有较少的线、面特征，故使用特征进行去畸变及配准精度可能性相对较低。本文以LOAM为基础构建了适合于森林条件的LiDAR SLAM算法，其工作流程如图2所示。

图2 森林LOAM工作流程图Fig.2 Workflow of LOAM for forest inventory

2.1 SLAM系统工作流程设计

针对森林中具有较少线、面特征的特点，对SLAM系统工作流程进行优化，以提高SLAM系统位姿估计鲁棒性并提高构建三维点云精度。相比于LOAM，本文构建系统进行了优化：

(1)在单帧点云特征提取时，将点云分为3类，即用于配准的下采样线特征、下采样面特征及其它点，在保存关键帧数据时，将用于配准的线、面特征存在内存中，其它点(占单帧点云总数约80%)则缓存至本地，相比于传统LOAM及LeGO-LOAM，该方案可构建密度更大、范围更大的三维点云。

(2)所使用线特征仅保留连续点形成的曲率较大点，遮挡形成线点被剔除(即仅保留两侧连续点与线特征点的深度差小于30 cm的特征)，防止立木切线被作为线特征参与运算。

(4)采用基于里程计的回环检测及基于表面的回环检测(Scan Context++[29])完成回环帧的判别，然后使用基于线、面特征的配准获取回环帧与当前帧约束；相比于ICP算法，基于线、面特征的特征关联方法更高效且可靠性更强。

(5)不仅为其构建了UI界面，通过进一步改造，使其摆脱ROS环境，可在Windows、Ubuntu等操作系统下编译并运行(图3)。

图3 森林LiDAR SLAM系统界面及其操作Fig.3 Forest LiDAR SLAM system interface and its operations

2.2 带有先验条件的SLAM优化算法设计

LOAM及其变体进行去畸变及配准优化中，将当前帧线/面特征与参考帧(或局部地图)特征进行关联后基于点到线、面的距离最小化完成位姿估计，点到线、面的距离可表达为

(3)

(4)

以上方程均可线性化为广义平差误差方程

Av=BX-l

(5)

其中

式中A——单位矩阵

v——点到线或面的偶然误差

B——线性化矩阵(即dL或dS对位姿线元素及角元素的导数)

l——给定初始位姿X0条件下点到线、面的距离负值

X——待求初始位姿的改正数

利用以上广义平差误差方程式，便可利用牛顿迭代或买夸特算法等估计待求改正数，从而逼近位姿真值，实现去畸变或位姿估计。

将激光雷达的先验条件等用于“去畸变”及“配准”算法，以提高位姿估计及建图精度。在“去畸变”优化中，线特征平差误差公式中参数可表达为

(6)

(7)

l=[-dL]

(8)

(9)

式中vdL——线特征到线的偶然误差

面特征平差误差公式中参数可表达为

(10)

(11)

l=[-dS]

(12)

(13)

式中vdS——面特征到面的偶然误差

(14)

转换为直角坐标系值。极坐标观测值的精度在激光雷达说明书中已经给出，或通过使用前检校可获取。若令

(15)

ΣX=MΣpMT

(16)

本文SLAM系统在“配准”时，将去畸变后线、面特征精度(协方差阵)代入到配准优化过程中，其中线特征配准平差误差公式参数可表达为

(17)

(18)

l=[-dL]

(19)

(20)

面特征配准平差误差公式参数可表达为

(21)

(22)

l=[-dS]

(23)

(24)

3 三维点云构建与立木检尺参数提取

3.1 研究区域概况

选择位于北京市近郊西山试验林场(39°58′N，116°11′E)为研究区域，该林场主要以人工林为主体。选择其中4块32 m×32 m的方形样地为研究对象。所选样地地面具有较少灌木且容易到达，不同径阶组比较均衡。表1总结了不同样地的基本属性。

表1 样地属性的描述统计Tab.1 Summary statistics of plot attributes

3.2 数据采集与处理

选择Velodyne VLP-32C型激光雷达为数据采集设备，该设备测距范围为200 m、典型场景下测量精度为±3 cm、垂直视场角为40°、水平视场角为360°、水平角分辨率最大可达0.1°(本文选择最高分辨率下进行样地扫描)。为方便手持并减少遮挡，为其加装了手持柄；为判断激光雷达竖直状态，在激光雷达顶端安装了万向水准器。改装后设备如图4所示。

图4 用于样地扫描的激光雷达系统Fig.4 LiDAR system for scanning forest plots1.水准仪 2.激光雷达 3.数据采集系统 4.手持柄 5.移动电源

为保证样地扫描点云完整性，并最大限度利用SLAM系统回环检测减少位姿漂移，以提高建图精度，设计了固定的样地扫描路径(图5)，即扫描起点为样地中心，沿西北方向开始扫描；到达西北角点后开始类航空摄影测量航线式扫描，扫描平行轨迹间距为6 m；完成类航线扫描后到达东南角，最后还需回到样地中心。修正扫描路径：①相邻平行轨迹间利用回环检测实现局部位姿漂移修正。②轨迹起始段及终止段轨迹与类航线相交，利用回环检测可实现局部位姿漂移修正。③轨迹终止点与起始点重合，利用回环检测实现全局扫描路径位姿漂移修正。

图5 样地扫描路径Fig.5 Plot scanning path

为保证激光雷达坐标系与其他参考坐标系的转换，具体样地数据采集过程为：①在样地中心及正北方向边界处放置激光反射标志(图6)。②手持激光雷达至样地中心，保持水平后开启数据采集。③沿规划路径完成样地扫描。

图6 激光反射标志Fig.6 Laser reflection mark

在完成森林样地扫描后，将数据导入本文构建LOAM系统完成样地三维点云构建；然后，使用LiDAR 360软件通过坐标变换、地面提取、DEM生成与高度归一化、胸高提取及胸高圆柱体拟合等，完成样地立木胸径及位置提取(图7)。

图7 森林样地点云后处理Fig.7 Post-processing of forest plot point cloud

3.3 研究方法

利用LiDAR SLAM系统及LOAM系统构建了森林样地三维点云。然后，将点云中提取立木胸径及立木位置与参考值对比，实现对SLAM系统生产点云精度的间接评估；通过对比SLAM系统与LOAM系统后处理数据统计结果，评价本文所构建SLAM系统。本文使用胸径尺测量胸径作为胸径参考值，全站仪(Leica TS60型)测量立木位置数据与胸径参考值联合计算立木位置为立木位置参考值。在使用全站仪对立木位置进行测量时，将全站仪架设于样地中心(即与布设于样地中心的激光反射标志十字中心对齐)，通过瞄准布设于样地正北方向的反射标志完成北向初始化。若由于遮挡，部分立木无法被观测，可通过“合并多站”方式完成所有立木位置观测。本文采用由立木位置误差统计获取的均值向量μ及二维协方差阵Σ对立木位置估计值精度进行评价，其定义式为

(25)

(26)

式中 (xi，yi)——立木位置估计值

(xir，yir)——参考值

n——立木总数

此外，协方差阵Σ的最大特征值σmax可表示为

(27)

式中 eigenvalues()——Σ所有特征值

该值为立木位置最大变异性方向的标准差，本文用其评价立木位置精度。

使用偏差(BIAS)、均方根误差(Root mean squared error, RMSE)、相对偏差(relBIAS)及相对均方根误差(relRMSE)对胸径估计值进行评估(其中RMSE也用于立木位置精度评估)。

4 试验结果分析

4.1 SLAM后处理数据定性评估

利用SLAM系统及LOAM系统分别对样地原始帧数据进行后处理，构建森林三维密集点云数据。从定性角度看：①2种SLAM解算获取的姿轨曲线并非完全重合，故2种解算方法获取森林点云具有一定区别(图8红色轨迹为森林SLAM解算姿轨曲线，蓝色轨迹为LOAM解算姿轨曲线；灰色点为样地点云)。②SLAM剔除了遮挡线特征，避免立木与视点切线被作为线特征参与运算，采用回环检测方法修正位姿图，建图漂移更小，不同关键帧拼接而成的胸高点云分层较LOAM系统小，故点云厚度较薄；而LOAM中将立木与视点切线特征作为线特征将产生误匹配，胸高圆柱体中心往往被填充(图9)。

图8 SLAM与LOAM后处理姿轨曲线定性对比Fig.8 Comparison of post-processing trajectory curves between forest SLAM and LOAM

图9 SLAM与LOAM后处理胸高点云数据定性对比Fig.9 Comparison of post-processing DBH point clouds between forest SLAM and LOAM

4.2 立木位置精度评估

由LiDAR SLAM系统及LOAM系统获取立木位置估计值，如图10所示，显然LiDAR SLAM估计值较LOAM偏差小。表2统计结果表明：①SLAM系统解算结果中，各样地x、y两轴估计值的平均误差为-0.029～0.027 m，即无明显偏差；而LOAM解算结果虽然总体接近无偏差，但各样地平均误差较大(-0.101～0.092 m)。②SLAM系统解算结果中，x、y两轴协方差值接近于0，说明两轴间估计值无明显相关性，LOAM解算结果有类似结论。③SLAM系统解算结果中，最大误差方向的协方差σmax(0.066～0.098 m)及RMSE(0.052～0.104 m)均明显小于LOAM计算结果(σmax为0.112～0.148 m，RMSE为0.114～0.144 m)，显然森林SLAM估计值变异性较小(图11)。

图10 立木位置图Fig.10 Estimated and reference stem positions

表2 立木位置估计值的精度Tab.2 Accuracies of stem position estimations

图11 所有样地立木位置误差Fig.11 Position errors of all trees in plots

4.3 胸径精度评估

由LiDAR SLAM系统及LOAM解算立木胸径估计值如图12所示，SLAM估计值较LOAM估计值更接近于两轴中线，显然SLAM胸径估计值可靠性更强。表3统计结果显示，SLAM解算各样地胸径估计值BIAS为-0.08～0.65 cm(relBIAS为-0.59%～2.42%)，与LOAM估计结果(BIAS为1.89～2.62 cm，relBIAS为10.09%～14.42%)相比具有较小偏差；RMSE为0.087～1.23 cm(relRMSE为3.61%～4.94%)，显然较LOAM估计值(RMSE为3.37～3.86 cm，relRMSE为14.48%～25.08%)精度更高。从图13可以看出，不同径阶SLAM立木胸径估计值均具有较小误差且变异性比较一致；而不同径阶LOAM胸径估计值虽变异性比较一致，但有约2 cm的估计偏差。

图12 胸径估计值散点图Fig.12 Scatter plot of estimated DBH

图13 不同径阶组胸径误差箱形图Fig.13 Errors of DBH observations under different DBH

表3 胸径估计值Tab.3 Accuracies of DBH estimations

5 结论

(1)针对森林环境中线、面特征少，无法精确建图的缺点，构建了LiDAR SLAM森林样地调查系统，以便仅利用单台多线激光雷达可精确完成森林样地调查。该SLAM算法利用二次去畸变、二次配准改善了位姿估计及点云地图的鲁棒性；该系统将激光雷达测量精度、位姿估计精度等先验信息融入去畸变及配准算法中，提高了位姿估计及点云地图精度。

(2)该系统在4块32 m×32 m方形样地中进行了测试，通过从该系统产生点云中提取的立木位置和胸径对系统精度间接评估。经与参考值进行对比表明，该系统所获取点云在立木位置及胸径估计时均能获取比LOAM算法精度更高的结果。显然，改进的LiDAR SLAM算法使单台多线激光雷达高精度完成森林样地调查成为可能。