基于RGB-D SLAM手机的森林样地调查系统研究

范永祥 冯仲科 陈盼盼 高 祥 申朝永

(1.北京林业大学精准林业北京市重点实验室， 北京 100083; 2.安徽农业大学理学院， 合肥 230036)

0 引言

森林具有多种环境功能和社会经济功能，与人类生存、发展息息相关。在环境方面，森林具有固碳、水土保持、调节气候、保持生物多样性等多方面的作用；在社会经济方面，森林为人类提供了食物、木材、能源等[1-2]。

森林资源信息能够为不同层面的决策者提供管理森林、计划造林等依据。森林清查的主要目的是收集不同决策者所需要的森林属性，通常调查属性包括生物量、蓄积量、树种等[3]。传统的地面清查方法是在调查区域合理布设样地，并对样地的森林属性进行清查，通过这些点估计来反演整个研究区域平均或总体的林业属性[4]。通常样地调查属性包括树种、胸径、树高，生物量、蓄积量等属性一般通过构建以树种、胸径、树高为自变量的模型进行估计[5]。除这些属性外，样地树木位置是解算林分空间结构参数如角尺度、混交度、大小比数等因子所必需的属性[6]。传统样地调查中，使用胸径尺或卡尺测量胸径，使用测高器测量树高，使用罗盘、皮尺或全站仪测量树木位置。这些方法不仅耗时耗力，而且无法克服观测者主观因素的干扰。Laser-relascope[7]、电子测树枪[8]等设备提供了同时测量树木胸径、树高、位置等一体化解决方案，但仍然无法克服安装复杂、观测难、主观因素干扰等问题。

激光雷达(Light detection and ranging, LiDAR)技术的出现及计算机运算能力的提高为样地清查提供了新的解决方案[9-11]。利用激光雷达对样地进行扫描并获取样地三维采样点云，然后从点云中客观地提取样地属性。地基激光扫描(Terrestrial laser scanning，TLS)作为一种地基激光雷达，已经被许多学者用于样地清查，并通过算法对森林属性进行提取和评估[5,12]。TLS通常以两种方式被使用：单次扫描和合并多次扫描。单次扫描即在样地中心建立站点并对样地进行扫描[13]。该方式数据处理简单，但在较密的森林中由于遮挡的原因经常发生漏测。多扫描即在样地中建立多个扫描站点，然后利用算法将多次的单扫描数据进行合并[14]。该方式虽然避免了遗漏问题，但合并单次扫描数据的难度较大。显然，TLS不适合在较大的样地中使用。移动激光扫描(Mobile laser scanning，MLS)的出现解决了TLS多次扫描数据难以合并的问题，从而可以在较大的样地中进行森林属性清查[15-17]。移动激光雷达通常为载有全球卫星导航系统(Global navigation satellite system，GNSS)、惯性测量单元(Inertial measurement unit，IMU)、激光雷达组合的运动平台。其中，GNSS+IMU的组合能够在移动激光雷达系统运动过程中对LiDAR的位姿进行测量，利用所测得的位姿可以对各个时刻的扫描数据进行合并。但在林下没有被GNSS覆盖的区域，MLS系统很难构建全局一致的点云。即时定位及构图(Simultaneous localization and mapping，SLAM) 技术的出现使得在林下无GNSS信号的区域定位成为可能。SLAM系统在运动过程中，利用单目相机、双目相机或RGB-D相机对周围环境进行观测，从而获取观测序列，然后使用该视觉观测序列对周围环境建图，并估计SLAM系统的位姿[18-19]。目前的研究中，将激光SLAM和双目SLAM作为MLS系统的辅助相对定位系统，从而使MLS系统在无GNSS信号样地中扫描样地时仍能获取全局一致的点云[20-22]，但MLS系统构造复杂、质量较大、价格昂贵，尤其对于轮式系统还需要较好的林下工作条件。

ToF (Time of flight) 相机[23-24]作为一种测量深度的消费级产品，具有体积小、价格便宜、功率低等优势，是激光雷达的一种替代品，且经常作为RGB-D SLAM的输入。随着智能手机运算能力的增强及SLAM算法的改进，RGB-D SLAM系统可以在手机上运行，从而使手机具有实时定位并构图的功能[25-26]。TOMATK等[27]使用这种手机扫描样地，并在线下从获取的点云中提取了样地树木位置、胸径。FAN等[28]使用这种手机实现了实时获取样地的单木属性——位置、胸径和树高，但并未评估样地属性。

本文利用带有RGB-D SLAM系统的手机构建森林样地调查系统，主要实现样地构建、每木检尺及林分/样地属性计算，并使用增强现实技术进行结果实时显示与交互。应用该系统在18块半径为7.5 m的圆形样地中进行测试，对估计结果进行评估。

1 理论与技术

1.1 SLAM理论

SLAM为以控制及观测为输入，在未知环境中未知位置构建周围三维地图的同时实时估计当前运动平台位姿的过程。即在SLAM过程中，运动平台位姿及所有路标点的位置在没有任何先验信息的条件下被实时估计。从概率分布的角度讲[29-30]，控制输入可以描述为运动模型P(xk|xk-1,uk)，其中，xk-1、xk为k-1、k时刻的位姿，uk为控制输入。

观测输入可以被描述为观测模型P(zk|xk,m)。其中，zk为k时刻的观测，m为所有路标点。

SLAM求解过程通常被视为具有马尔可夫性，即未来的位姿仅与前一时刻的位姿有关系而与过去的位姿没有关系，从而运动模型与观测模型可以通过贝叶斯过滤(Bayes filter)来求解当前位姿的后验分布，先验位姿的估计被描述为预测更新

(1)

式中Z0:k-1——0到k-1时刻的观测

U0:k——0到k时刻的控制输入

x0——初始时刻的位姿

利用先验分布及观测模型求解后验位姿估计被描述为观测更新

(2)

贝叶斯过滤仅从概率角度解决了SLAM问题，在实际使用时需要给出运动模型及观测模型的具体形式，目前主流的解决方案包括扩展卡尔曼滤波(Extended Kalman filter，EKF)及其变体[31]、非线性优化(Sparse non-linear optimization method)[30]及粒子滤波(Particle filter)[31]。

1.2 RGB-D相机技术

RGB相机是一种用于获取周围环境纹理信息的工具，常用于SLAM系统的观测输入，但由于其仅具有测角而不具有深度测量的功能，所以无法为SLAM系统所估计的位姿及地图提供尺度。深度相机因具有测距功能，弥补了RGB相机不能获取深度的问题。RGB-D相机作为RGB相机与深度相机硬连接组成体，可以同时获取周围环境的纹理信息及深度信息。

目前,深度相机进行深度估计主要基于飞行时间(Time of flight)原理，这种相机被称为ToF相机。飞行时间一般通过使用脉冲信号或调制连续波来获取，但目前主要使用基于调制连续波的技术。这种技术通过使用发光二极管(Light-emitting diodes, LED)向目标发射近红外光(Near-infrared light, NIR)，当光被反射回来时利用传感器接收并获取当前相位差及信号强度，用于计算测量目标深度及其可靠性。ToF相机测距原理基本与激光雷达一致，但激光雷达采用逐点扫描的方式对周围目标深度进行采样，即在测量深度时需要记录角度才能获取所采样点的坐标；而ToF相机使用普通相机中心投影、通过像素坐标的方式测量角度。

将ToF相机与RGB相机硬连接后的组合体标定其内方位元素及相对外方位元素后便可以成为真正意义的RGB-D相机，然后便可将深度相机获取的像素深度映射到RGB相机像素对应的深度

(3)

其中

ud=[udvd1]T

(4)

(5)

(6)

式中dd——深度相机像素点坐标ud对应的深度

Kc、Kd——RGB相机与深度相机的内方位元素矩阵

fx、fy——x、y方向以像素为单位的焦距

cx、cy——投影中心投影在像素平面坐标系的坐标

Rc_d、tc_d——深度相机坐标系相对于RGB相机坐标系的旋转矩阵及平移向量

dc——RGB相机像素点坐标uc对应的深度

ud、uc——深度相机像素点坐标、RGB相机像素点坐标

1.3 Google Tango技术

虽然近年来SLAM算法得到改进，芯片运算能力得到提高，要使移动设备能够实时估计位姿仍然具有难度。Google Tango通过使用特殊硬件使智能手机能够运行SLAM系统。如图1所示，Google Tango以运动跟踪相机(RGB相机)及ToF相机(深度相机)作为观测输入，9轴的加速度传感器、重力传感器、罗盘传感器估计控制输入，并使用专用的低功率计算机视觉处理器加快数据处理过程。Google Tango手机使用基于外观的实时构图(Real-time appearance-based mapping，RTAB-Map)[32-33]SLAM系统，如图2所示，RTAB-Map系统分为前端和后端两部分。前端为视觉里程计，该部分使用RGB-D及IMU获取的数据逐帧估计手机的相对位姿。随着估计过程的进行，估计误差累积使位姿漂移逐渐加重，后端的目的正是为了减少这种误差，后端由局部光束法平差(Local bundle adjustment)模块与回环检测模块组成， 局部平差模块用于对最近的关键帧位姿进行平差修正；回环检测通过检测访问过的区域使用位姿图优化(Pose graph optimization)对各时刻的位姿进行纠正。

图1 Google Tango手机(Lenovo Phab 2 Pro)的结构Fig.1 Structure of Google Tango phone (Lenovo Phab 2 Pro)

图2 RTAB-Map系统位姿估计、地图构建过程Fig.2 Process of estimating pose and mapping in RTAB-Map system

2 样地调查系统设计与原理

该样地调查系统主要依托于SLAM系统。如图3所示，SLAM系统以RGB-D相机及IMU为输入并实时估计获取的实时RGB图像、点云数据、位姿等信息；样地调查系统使用位姿及深度信息来估计单木参数，RGB图像及位姿则通过OpenGL构建三维虚拟场景，通过将SLAM坐标系统与OpenGL坐标系统对齐从而通过手机屏幕为使用者提供增强现实(Augmented reality, AR)场景，并使用户能通过屏幕与AR场景进行交互。

图3 样地调查系统基本结构Fig.3 Structure of forest inventory system

样地调查系统主要包含定义样地坐标系、构建全局一致的稀疏地图、每木检尺及参数计算4部分，如图4所示。定义的样地坐标系用于描述样地中每一棵树木的位置；构建全局一致性稀疏地图使每木检尺时所获取的手机位姿通过回环检测减少漂移，进而约束树木位置的估计误差；每木检尺过程观测样地中的所有树木；参数计算过程计算样地所代表区域的林分参数等。

图4 样地调查系统工作流程图Fig.4 Workflow of forest inventory system

SLAM系统会定义一个以系统启动时设备位置为原点、x轴正向水平向右、z轴正向垂直向上、y轴正向指向深度方向的右手坐标系，即初始化坐标系。SLAM系统的位姿将以该坐标系为参考坐标系。样地坐标系是以样地中心为原点、x轴水平向东、y轴正向指向正北、z轴垂直向上的右手坐标系,可以更好地描述样地中树木位置。因此，需要获取两个坐标系的相对关系以实现将以初始化坐标系为参考系的树木位置转换为以样地坐标系为参考系。两坐标系之间的相对关系由平移和旋转共6个自由度来确定，仅需获取2个点在不同坐标系中的坐标值便可以获取两坐标系的变换矩阵。本文通过采集样地中心点xi_center以及正北方向一点xi_north来获取变换矩阵,其旋转矩阵Ri_p及平移矩阵ti_p分别为

(7)

ti_p=xi_center

(8)

(9)

为了自动获取样地中心及正北方向一点，本文设计了棋盘如图5所示，用其中3×3圆形格网指示样地中心、4×4棋盘格中心指示样地正北方向点，然后通过使用OpenCV库棋盘格检测的方法检测两棋盘格并获取圆形格网及棋盘格中心坐标，从而获取xi_center及xi_north，并经过转换获取样地坐标系(图6a)。

图5 确定样地坐标系所使用的棋盘格Fig.5 Chessboard used for building plot coordinate system

SLAM系统通过回环检测及位姿图优化可以减少前端视觉里程计位姿估计产生的漂移。若在观测前获取了样地全局一致的地图，则在每木检尺时手机位姿漂移将通过后端得到优化，所以本文在每木检尺之前进行了样地扫描从而构建样地稀疏地图的步骤。全局一致性地图的准确性与扫描轨迹有着密切联系，扫描过程适当的应用回环检测则可获取准确的地图，本文采用螺旋形路径进行扫描以尽可能多的构建回环检测(图7)，图中绿色三角形表示样地中心，也是扫描起点；红色方形为扫描终点，该点接近于样地边界。为了确保能够获取完整的样地地图且不至于浪费时间扫描样地外的地面， 或者在每木检尺时避免采集样地外树木的信息，确定样地边界极为重要。本文使用增强现实技术在OpenGL坐标系内构建了样地边界，当观测者接近样地边界时，手机屏幕将会显示样地边界的位置(图6b)。

图6 样地观测过程中样地调查系统状态Fig.6 Different statuses of plot survey system during observation

每木检尺过程需要对样地中所有立木的胸径、树高和树种等进行观测。为确定胸径所在位置，首先需要确定地径处一点以确定胸高位置(图6c、6d)，然后利用胸高处的点云过滤、拟合可获取立木胸径及位置信息[28]，经过映射到OpenGL坐标系便可以将结果展示在屏幕上(图6e)；在完成胸径、位置估计后，若需要观测树高则移动到能够观测到树梢的位置，在屏幕上点击树梢，便可计算树高 (图6f)。

图7 样地扫描路径Fig.7 Scanning trajectory for building plot map

每木检尺过程获取了立木的胸径、树高、立木位置、树种信息(图6g、6h)，所有样地、样地所代表的林分等信息便可由此计算(图6i)。林分平均胸径为

(10)

式中di——样地中第i棵树的胸径，cm

N——样地中立木总数

林分平均树高为

(11)

式中hi——样地中第i棵树的树高,m

第i棵立木材积为

(12)

式中fε——实验形数

每公顷蓄积量(单位：m3/hm2)为

(13)

式中S——样地面积，m2

每公顷横断面积(单位：m2/hm2)为

(14)

林分株树密度(单位：株/hm2)为

(15)

若样地中立木的位置平均值μ及协方差Σ分别为

(16)

(17)

其中

xi=(xi,yi,zi)T

式中xi——样地中第i棵树在样地坐标系中的位置坐标

则样地近似平面的法向量为协方差矩阵Σ最小特征值所对应的特征向量nmin_eigenvalue，则样地的坡度及坡向分别为

(18)

(19)

其中

(20)

式中Iz——样地坐标系中z轴方向的单位向量

Iy——样地坐标系中y轴方向的单位向量

E——将nmin_eigenvalue投影到Oxy平面的投影矩阵

3 实验区域与方法

3.1 实验区域概况

实验区域为北京市西山实验林场，该林场位于北京市近郊小西山(39°58′N，116°11′E)。选择18块半径为7.5 m的圆形样地为研究对象。为了使实验具有代表性，所选的样地具有不同的坡向、坡度，林分具有不同密度、树种。所选样地地面具有较少灌木且容易到达。表1为不同样地立木的基本属性。

表1 样地属性Tab.1 Summary statistics of plot attributes

3.2 实验方法

为了验证设计的系统在样地调查中的精度，使用胸径尺测量胸径作为胸径参考值，使用全站仪(Leica Flexline TS06plus)测量树高作为树高参考值，使用全站仪测量立木位置作为立木位置的参考值。样地所代表区域林分参数的参考值均基于胸径、树高及立木位置参考值进行计算。为了给出所设计系统的估计精度，本文使用偏差(BIAS)、均方根误差(Root mean squared error, RMSE)、相对偏差(relative BIAS)及相对均方根误差(relative RMSE)对各测量值进行评估，计算公式如下

(21)

(22)

(23)

(24)

式中xi——测量值

xir——测量值xi对应的参考值

n——测量值总数

BIAS——估计值偏差

Brel——估计值相对偏差

RMSE——估计值均方根误差

Rrel——估计值相对均方根误差

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

4 实验结果

所设计系统通过样地调查估计林分因子及样地因子的关键在于测量样地中立木胸径、树高及位置的精度。如图8所示，胸径估计值是无偏的且随着胸径增加估计值的误差离散度不变。由表2可以看出，胸径估计值具有较小的偏差(0.32 cm)及较小的离散度(1.27 cm)。图9说明了树高估计值是无偏的(8.0 cm)，而估计值的离散度随树高增加而变大，但由表2可知总体仍然具有较小的离散度(95 cm)。图10描述了树木在水平面上的位置误差，总体而言立木位置是无偏的(x方向偏差为3.0 cm；y方向为-1.7 cm)，且离散度较小(x方向为3.8 cm；y轴方向为3.0 cm)，由表2可知，高程方向的偏差(0.4 cm)及离散度(1.1 cm)相较于x、y方向较小。

表2 单木属性精度Tab.2 Estimations of tree-based properties

图9 树高估计值散点图Fig.9 Scatter plot of estimated tree height values

图10 立木平面位置估计误差Fig.10 Position errors of all trees in plots

表3为18块样地的林分参数，所有的林分参数均与参考值相接近。仅样地3株树密度与参考值不一致。如表4所示，总体而言所有林分估计值具有较小偏差及较小的离散度。表5比较了样地坡度及坡向测量值与估计值，结果表明坡度的估计值与参考值接近；当坡度较小时，由于SLAM系统位姿估计误差导致坡向估计值与参考值差异较大。总体而言，坡度及坡向的测量值仍然可靠。

表3 林分属性估计结果Tab.3 Estimation results of forest-based properties

表4 林分/样地属性估计值的精度Tab.4 Accuracies of forest-based and sample-based properties

5 结论

(1)基于RGB-D SLAM手持式设备构建了森林样地调查系统，该系统使用SLAM系统估计了设备的相对位姿，利用ToF相机获取周围三维采样，依此设计了计算立木胸径、树高及位置的算法，从而形成了构建样地、每木检尺及林分/样地参数估计的调查系统。基于Android平台显示系统设计了增强现实功能，用于实时显示调查结果，为观测者提供检查及交互功能。

表5 样地属性估计结果Tab.5 Estimation results of plot-based properties (°)

(2)该系统在18块半径为7.5 m的圆形样地中进行了测试。与参考值对比表明，该系统在林分及样地属性估计时均能获取较精确的结果。