基于视觉里程计的森林样地调查系统研究

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

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

0 引言

森林占全球陆地生态系统的30%，为全世界50%以上的动植物提供了栖息地，是生物遗传多样性最为丰富的生态系统[1-5]，在气候调节、碳循环和水土保持等方面发挥着重要作用[6-7]。森林清查是自然资源监测和管理的重要组成部分，是了解森林资源现状及其动态变化的主要方式[8]。森林清查依不同层次的森林管理者及政府决策者的需求不同而不同，为森林蓄积量、生长量、生物量及碳储量的分析和评估提供系统的森林资源信息[9]。

森林样地调查是森林资源清查的重要方式之一,不仅可以直接用于评估森林资源信息，也是空基等数据构建模型的参考数据源。森林样地通常为森林中具有代表性的小面积区域，通过对样地中的树木进行每木检尺获取属性信息并加以总结，便可获取森林总体的评价属性。森林样地调查的单木属性主要包括树种、位置、胸径、树高等，而生物量、蓄积等属性可通过所构建模型进行估计[10-11]。显然，样地调查的精度及效率直接影响森林清查的可靠性及效率[12]。在传统森林样地调查中，单木属性信息的获取主要依赖于简单的测量工具，如胸径尺、测高器、测树枪等。这些方法耗时、费力，且存在较大主观因素[13-15]。

随着遥感技术的发展，利用传感器获取样地点云，然后从该点云中提取立木参数的方式被许多研究者关注[16]。其中，主要使用的传感器包括激光雷达(Light detection and ranging， LiDAR)和摄影测量(Photogrammetry)相机。激光雷达通过主动测量目标点与激光雷达距离的方式获取深度信息[17]。由于遮挡等原因导致该方法测得的树木点云密度过小，甚至无法进行拟合以获取树木胸径及位姿[18]。摄影测量是一种被动对周围纹理进行观测的工具，通过匹配、位姿估计、稠密化点云等步骤获取森林三维点云。但由于森林结构复杂，不一定能够成功获取同名点，而获取稠密的树木点云，当然也无法获取较高精度的单木属性[19-22]。

利用摄影测量很难保证所获稠密点云合格并可以提取所有立木的单木属性[23-28]。视觉里程计是一种仅使用相机便可以基于摄影测量理论实时估计相机位姿的一种技术。显然，这种不需要全球导航卫星系统(Global navigation satellite system，GNSS)信号覆盖条件的技术适合在林下工作。本文设计一种以视觉里程计估计森林中获取的连续图像序列的位姿、利用非点云的方式获取立木胸径及位置的方法。利用该方法构建图像序列处理系统，并将其用于样地中获取图像序列处理，以获取样地中立木胸径及位置。

1 视觉里程计理论

视觉里程计(Visual odometry，VO)为以搭载在运动平台上的视觉传感器作为输入，利用传感器获取的时间序列帧数据，进行实时估计搭载平台位姿的技术[29]。VO系统所搭载的视觉传感器主要包括单目相机(Monocular camera)、双目相机(Stereo camera)和深度相机(RGB-D camera)3类[30]。视觉里程计工作所需的传感器较为简单，而且在不需要GNSS覆盖的条件下便可进行相对位姿估计，可以用于机器人室内相对导航及宇宙探索中，也适合于林下无GNSS信号条件下估计观测设备相对位姿的情况[31]。

VO算法主要需要解决的问题是利用不同视觉传感器获取的具有时间相关性的图像序列进行实时估计搭载平台的位姿。若设视觉传感器在N个不同时刻获取的帧序列为I0:N={I0,I1，…,In}(图1)，视觉里程计通过相邻帧对同一环境进行重叠观测的特性从而估计两相邻帧k-1与k之间位姿的相对变换，该刚性变换可以表示为

(1)

式中Rk-1,k——从k-1时刻到k时刻的旋转矩阵

tk-1,k——从k-1时刻到k时刻的平移向量

所有连续帧的相对变换集合为T1:n={T0,1,T1,2,…,Tn-1,n}，若设不同时刻的位姿序列为C0:n={C0,C1,…,Cn}，则k时刻的位姿为

Ck=Ck-1Tk-1,k

(2)

且通常设k=0时刻的位姿C0作旋转矩阵为单位阵、平移向量为零向量的假设。基于此便可以估计不同时刻的位姿。

图1 视觉里程计的本质Fig.1 Essence of visual odometer

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

2 森林样地数据处理系统设计

森林样地数据处理系统主要包括视觉里程计及森林调查系统两部分(图2)。视觉里程计部分以外业所采集的图像序列为输入估计各帧图像的位姿；森林调查系统需要手动标记树木位置及胸径，然后基于同名点自动匹配构建样地坐标系并估计立木位置和胸径。

2.1 视觉里程计

通常使用的视觉里程计为基于特征的里程计，其利用不同帧中出现多个相同特征构建不同帧之间的位姿转换关系。如图3所示，为本文所构建的基于特征点的视觉里程计估计位姿的基本流程图，其利用有序的图像序列，通过特征提取、特征匹配、运动估计和局部光束法平差4个步骤为获取当前帧的位姿：① 特征点提取部分用于提取图像中对应于地面的点，好的特征点应具有位置精确(位置及尺度)、可重复性、计算高效及强鲁棒性等特征，本文选择提取速度较快的FAST(Feature from accelerated segment test)方法，提取角点类型的特征点。②特征匹配部分需要建立不同帧中特征点之间的对应关系，并使用BRIF(Binary robust independent elementary feature)算子描述该特征点的特点且该算子之间的距离可以描述两个特征点的相似性，便可依次获取不同图像、相同场景特征点的对应关系。③当获取到匹配点时，即可利用P2P或EPNP估计当前帧的位姿。④利用局部光束法平差从统计优化的角度精确化所估计的位姿。

图3 视觉里程计工作原理Fig.3 Operating principle of visual odometer

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

图4 基于全局地图的视觉里程计Fig.4 Visual odometer based on global map

为了减少位姿估计过程中的漂移，本文视觉里程计采用全局稀疏特征点地图的方式进行研究(图4)。特征点地图中特征包含描述子、位置及其标准差信息。当获取到第1帧图像的特征点时，假设该帧所在位置为世界坐标原点且无旋转，故直接将所有特征点添加到特征点地图中并赋予适当的位置及标准差。当获取第2帧数据时与地图中的特征点进行匹配，基于匹配的点利用P2P算法、局部光束法平差估计该帧的位姿，获取到位姿后基于三角测量重新估计特征点地图中特征点的位置及其标准差；除此之外应将新获取的特征点添加到特征点地图。此后获取的帧数据则使用EPNP及局部光束法平差估计该帧图像的位姿，在更新特征点部分与第2帧相同。一个典型的位姿估计结果如图5所示。

图5 位姿估计结果示例Fig.5 Example of pose estimation results

2.2 立木位置与胸径估计系统构建

利用样地中拍摄的有序图像数据估计样地中立木的位置及胸径，通常所需要的立木位置使用以x轴水平向东、y轴正方向指向正北、z轴垂直向上的右手坐标系(样地坐标系)进行描述。然而，视觉里程计的全局坐标系是以第1帧图像的位姿构建的初始化坐标系。另外，视觉里程计假设第2帧图像与初始化坐标系原点的距离为单位距离。显然视觉里程计初始化坐标系描述下的位姿尺度不正确，且与目标坐标系——样地坐标系不一致。显然，在估计样地中立木位置前需要纠正其尺度，并构建从初始化坐标系到样地坐标系的转换，从而使所有帧数据能够在样地坐标系中描述，从而提取基于样地坐标系下的立木位置。然后，不同于传统基于MVS(Muti-view stereo)构建稠密点云并利用点云提取立木信息的方法，本文试图利用在图像上直接标注立木位置及胸径的方法估计待估单木信息。

2.2.1样地坐标系构建

为构建样地坐标系，在扫描样地前首先在样地中放置如图6所示的特殊棋盘格，并使棋盘格上对应于样地中心的点与样地中心重合，箭头方向指示正北方向。在内业处理时，按照软件提示选择相关帧棋盘格对应的样地中心点及正北一点(图7a)，然后使用极线搜索和块匹配技术[32]搜索其他相关帧对应于这两点的像素坐标，并通过三角测量的方法估计这两点在初始化坐标系的坐标值，然后便可依据棋盘格对应这两点的实际距离修正所有帧的实际位姿。尺度计算公式为

(3)

式中L——棋盘格上两点的实际距离

xnorth、xcenter——棋盘格上正北一点及中心点在初始化坐标系中的坐标值

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

在修正尺度后便可获取正确尺度的样地中心及正北一点，为确定样地坐标系与初始化坐标系之间的变换关系还需要获取竖直方向的一个向量，本文通过假设样地中立木树干竖直向上估计该向量。如图7b所示，通过构建的森林样地调查系统在图像中选择立木地径及树干上一点，并使用极线搜索、块匹配技术搜索自动估计相应的在初始化坐标系中的坐标，则竖直方向的单位向量为

(4)

式中xbottom_i、xtop_i——第i棵被标记立木的地径处一点及树干上一点的坐标

N——被标记立木总数

由此，便可以估计样地坐标系x轴单位向量在初始化坐标系中为

(5)

y轴的单位向量为

Iy=IzIx

(6)

由此可以获取由样地坐标系到初始化坐标系的旋转矩阵为

Rinit_plot=[IxIyIz]

(7)

及平移向量

tinit_plot=xcenter

(8)

可以将所有图像位姿以样地坐标系进行表示。

2.2.2立木位置及胸径估计

图7 森林调查系统操作界面Fig.7 Forest survey system operation interfaces

(9)

式中xplot——胸高处一点在样地坐标系中的坐标值，其可由地径处点坐标z轴值加1.3 m获取

Rcamplot、tcamplot——图像由样地坐标系变换到相机坐标系的旋转矩阵及平移向量，即图像在样地坐标系中的位姿的逆变换

K——相机内方位元素矩阵

zcam——xplot在相机坐标系下的z轴坐标值

在获取到胸高处一点，基于极线搜索和块匹配精确化胸径处点，确保该点落在树干上。然后，需要在森林调查系统软件中点击树木胸高处左、右侧边缘以完成树木位置及胸径的估计(图7d)。在点击立木左右边缘后，一个以该图像的位姿原点为原点、z轴水平指向胸径处点、y轴竖直向上的右手辅助相机坐标系被构建以便于参数估计，从而求解问题可以被描述到oAC-xAC-zAC平面内(图8)。则立木的平面位置及胸径计算式为

(10)

(11)

其中xC=[xCzC]TxL=[xLzL]T

xR=[xRzR]T

图8 树木位置及胸径估计问题Fig.8 Tree position and breast diameter estimation problem

式中xC、xL、xR——树木位置在oAC-xAC-zAC平面的坐标

d——胸径

xL、xR——树干边缘处点的平面坐标

xcorr_bh——xL和xR的中心点坐标

其由胸径处点xp=[xpzp]T及左右侧边缘的像素获取得

(12)

式中Rac_c——从相机坐标系到辅助坐标系的旋转矩阵

在获取到立木位置的平面位置xC后设置其高度为胸径处点的高，便可换算到样地坐标系从而获取到样地坐标系描述下的立木位置。

3 实验材料与方法

3.1 研究区域概况

以北京市西山实验林场(39°58′N，116°11′E)为研究区域。该林场海拔为300～400 m，平均降水量为630 mm。研究区主要包括人工纯林和针阔混交林。选择了12块半径为7.5 m的圆形样地为研究对象，被选择样地地面灌草较少，以尽量减少环境对相机视野的阻挡。样地的基本情况如表1所示。

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

3.2 实验方法

采用大疆灵眸Osmo型口袋云台相机，其成像传感器为1/2.3”CMOS，该相机具有质量小、结构紧凑、智能、便携等优势，另配置有较小的三轴机械增稳云台，抖动抑制精度高达±0.005°。通过Matlab 2014a软件，利用张正友标定法对相机进行标定，所获取相机内参及畸变系数如表2所示。

为了验证森林样地内业处理系统所估计单木参数的精度，使用胸径尺测量胸径作为胸径参考值，立木位置使用全站仪测量值联合胸径计算获取的立木位置作为树木位置参考值。本文使用偏差(BIAS)、均方根误差(Root mean squared error, RMSE)、相对偏差(relBIAS)及相对均方根误差(relRMSE)对各观测值进行评估，其计算公式为

表2 相机内参标定数据Tab.2 Experimental result for camera calibration of intrinsic parameters

(13)

(14)

(15)

(16)

式中xi——测量值

xir——相对于测量值xi的参考值

n——测量值的总数

BIAS——估计值偏差

Brel——估计值相对偏差

RMSE——估计值均方根误差

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

4 结果与分析

4.1 立木位置估计

立木位置估计值在x、y轴方向分别为0.04、-0.03 m，且最大偏差不超过0.10 m(表3)。虽然样地6在x轴方向的RMSE达到0.40 m，但总体而言估计值的RMSE(0.21、0.17 m)较小。图9为所有样地中立木位置的误差，由图9可知，虽然单个样地中似乎立木位置存在系统误差，但总体而言其误差趋于0 m，且最大误差小于0.7 m。

4.2 胸径估计

经与参考值对比，胸径估计值具有较小偏差(0.09 cm, 0.51%)及较小的RMSE(0.88 cm，5.03%)(表4)。其中，最大偏差及RMSE分别为-0.95 cm及1.16 cm。由图10可见，胸径估计值与参考值相接近，且无偏差过大点；且不同径阶估计值也无明显偏差且具有相近的变异性(图11)。

5 讨论

首先利用视觉里程计基本原理恢复了在森林样地中获取图像序列的位姿，然后利用构建的森林样地调查系统通过构建样地坐标系、三角测量等过程估计了样地中立木胸径及位置。测试结果表明：立木位置估计值在x轴和y轴方向的BIAS分别为0.04、-0.03 m，RMSE分别为0.21、0.17 m；立木胸径估计值的BIAS及RMSE分别为0.09 cm(0.51%)和0.88 cm(5.03%)。

表3 立木位置估计值的精度Tab.3 Accuracies of tree position m

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

样地号RMSE/cmrelRMSE/%BIAS/cmrelBIAS/%10.914.900.442.3320.673.130.080.3530.593.880.130.8840.472.500.160.8551.166.820.633.6761.086.680.875.3570.854.01-0.25-1.1880.895.190.362.0791.047.020.392.65100.844.890.241.42111.046.42-0.95-5.85120.844.92-0.71-4.17平均0.885.030.090.51

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

图11 胸径估计值箱线图Fig.11 Errors of DBH observations for different DBH values

目前对立木位置的研究相对较多，但较少明确给出精度。TANG等[32]利用移动激光雷达对不同的样地进行了扫描，其运动轨迹的RMSE为0.32、0.29 m，但未评估树木的精度。FAN等[10]利用RGB-D SLAM手机对立木位置进行了实时估计，其x轴及y轴方向均方根误差(RMSE)均为0.12 m。LIU等[33]设计了一种实时动态多功能CCD连续摄影的地面观测仪器，通过点云提取立木位置，其偏差最大为0.19 m，最小为0.137 m，均方根误差(RMSE)最大为0.201 m，最小为0.162 m。显然，本文所构建调查系统所获取立木位置精度与以上方法所获取结果相似。

目前对立木胸径获取方法的研究较多。LIANG等[34]通过地面激光雷达获得点云数据并利用圆柱体拟合的方法来提取胸径，其胸径估计值的BIAS为-0.18～0.76 cm，RMSE为0.74～2.41 cm；范永祥等[35]基于RGB-D SLAM手机所构建的森林样地调查系统对立木胸径进行了评估，结果表明：胸径估计值的BIAS及RMSE分别为0.36 cm(1.93%)及0.69 cm(3.89%)；QIU等[36]利用连续摄影测量原理及图像识别技术研发了PDA(Personal digital assistant)摄影测树仪器，该设备估计胸径最大偏差为2.1 cm，绝大多数样地测量偏差小于1.5 cm。本文所构建系统所获取胸径精度与这些设备也相似。

6 结论

(1)设计了以视觉里程计估计、从森林样地中获取的有序图像序列的位姿；并构建了森林样地数据处理系统，以实现尺度修正、建立样地坐标系、立木位置及胸径标记、获取样地中立木位置及胸径。

(2)选取12块半径为7.5 m的圆形样地为研究对象，利用相机获取样地中连续图像序列。经所设计系统处理，所获取胸径及立木位置均具有较高精度。